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A  generalized,  finite  volume-based  SOFC  model  is  presented  which  includes  charge,  mass  and  heat  trans¬ 
port  as  well  as  a  Langmuir-Hinshelwood  type  applied  kinetic  model  for  steam  reforming  reactions.  The 
model  development  aimed  at  fast  applicability  to  various  cell  geometries,  short  calculation  times  to  allow 
for  system  analysis  calculations  and  high  fuel  flexibility.  In  the  first  part  of  the  paper,  the  model  approach 
and  assumptions  are  presented  in  detail.  In  the  second  part,  the  generalized  model  is  applied  to  a  pla¬ 
nar,  the  standard  tubular  and  the  triangular  tube  cell  geometry.  The  validation  against  experimental  and 
benchmark  test  data  is  stressed  to  assure  comparability  of  the  model  results  for  the  three  investigated 
cell  designs.  In  the  last  part,  the  three  cell  designs  are  compared,  highlighting  the  differences  with  respect 
to  internal  heat-  and  mass  transfer  and  the  impact  on  the  electrochemical  performance.  It  is  shown  that 
the  performance  of  the  triangular  tube  cell  is  almost  double  that  of  the  standard  tubular  cell  designs.  The 
planar  cell  is  outperformed  by  the  triangular  tube  cell  as  well. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cells  (SOFC)  feature  high  efficiency  in  convert¬ 
ing  gaseous  fuels  to  electricity  at  scales  ranging  from  kilowatts  to 
megawatts.  Hydrogen  and  carbon  monoxide  are  fuels  for  SOFCs 
and  the  operational  temperatures  departing  between  700  °C  and 
1000  °C  enable  the  indirect  conversion  of  hydrocarbons  via  steam 
reforming  reactions.  Recently,  evidence  for  direct  conversion  of 
hydrocarbons  was  found  [1].  For  this  reason  SOFCs  have  widely 
been  considered  as  the  most  fuel-flexible  fuel  cell  type.  Besides, 
the  high  exhaust  temperatures  predestine  SOFCs  for  the  combina¬ 
tion  with  gas  turbines  (GT).  Such  SOFC/GT  systems  are  expected  to 
reach  electrical  efficiencies  up  to  65%  [2]. 

SOFC  development  is  still  in  progress  and  a  number  of  cell 
designs  with  specific  advantages  and  drawbacks  are  under  inves¬ 
tigation.  The  most  important  cell  designs  are  the  planar  and  the 
tubular  designs.  Planar  SOFCs  feature  short  current  paths  result¬ 
ing  in  low  electrical  resistance  and  high  power  densities.  However, 
planar  cells  suffer  from  thermal  stress  resulting  from  different 
thermal  expansion  coefficients  of  the  employed  materials  for  the 
anode-electrolyte-cathode  (AEC)  assemblies,  the  interconnector 
plates  (IC)  and  the  sealings.  Tubular  cells  have  considerably  less 
problems  with  thermal  stress  and  are  therefore  very  robust.  Long- 
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duration  tests  with  tubular  cells  up  to  69,000  h  were  reported  [3]. 
The  seal-less  tubular  cell  design  involves  however  long  current 
paths  and  higher  electrical  resistance  than  planar  cells.  The  current 
path  length  further  limits  the  size  of  tubular  cells  which  compli¬ 
cates  a  cost-effective  production.  The  triangular  tube  cell,  in  the 
following  referred  to  as  Delta8  cell,  combines  the  advantages  of  the 
seal-less  design  with  shortened  current  paths  and  larger  cell  sizes 
which  result  in  lower  production  cost. 

From  the  mentioned  differences  of  these  three  cell  designs,  the 
question  arises  what  the  impact  of  the  specific  cell  design  character¬ 
istics  on  the  overall  system  is.  This  point  is  even  more  important, 
when  non-standard  fuels  such  as  producer  gases  containing  tars 
and  hydrocarbons  shall  be  applied.  Purely  experimental  investiga¬ 
tions  are  hardly  applicable  to  solve  this  issue  due  to  the  high  costs 
involved.  Therefore,  numerical  modeling  of  SOFCs  has  drawn  a  lot 
of  interest  in  the  past  years. 

Usually  simplified  models  have  been  employed  for  systems  anal¬ 
ysis  calculations,  concentrating  on  the  overall  cell  performance, 
e.g.  [4-8].  This  approach  is  valid  at  standard  operating  conditions; 
however  it  is  questionable  for  the  prediction  at  non-standard  oper¬ 
ating  conditions  such  as  e.g.  operation  with  other  fuels  or  electrical 
loads.  For  such  investigations,  more  complex  models  are  required 
which  consider  the  internal  charge,  heat  and  mass  transfer  pro¬ 
cesses  locally  resolved  instead  of  averaged  over  the  entire  cell. 

Typically  computational  fluid  dynamics  (CFD)  codes  have  been 
employed  for  cell  design  optimization.  The  calculation  times  of 
25-50  min  [9],  however,  hindering  an  application  of  CFD-based 
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Nomenclature 

A 

area  (m2) 

AEC 

anode-electrolyte-cathode  assembly 

cp 

heat  capacity  (J  mol-1  K-1 ) 

CV 

control  volume 

d 

diameter  (m) 

D 

diffusion  coefficient  (m2  s-1 ) 

Dk 

Knudsen  diffusion  coefficient  (m2  s-1 ) 

Dm 

molecular  diffusion  coefficient  (m2  s-1 ) 

E 

voltage  (V) 

^act 

activation  energy  for  exchange  current  density  cal¬ 
culation  (J  mol-1 ) 

DQh2 

hydrogen  equivalent  molar  flow  (mol  s-1 ) 

F 

Faraday  constant  (C  mol-1 ) 

h 

height  of  gas  channel,  bipolar  plate,  etc.  (m) 

I 

current  density  (Am-2) 

h 

exchange  current  density  (Am-2) 

Kp 

equilibrium  constant 

l 

length  (m) 

k 

length  of  planar  cell  (m) 

b8 

length  of  Delta8  cell  (m) 

h 

length  of  tubular  cell  (m) 

L 

length  of  cell-subunit  (m) 

LHV 

lower  heating  value  (J  mol-1 ) 

M 

molecular  weight  (kg  mol-1 ) 

n 

number  of  exchanged  electrons  per  electrochem. 
reaction 

nc  h 

number  of  channels  of  planar  cell 

ndot 

molar  flow  (mol  s-1 ) 

Hre 

number  of  repeating  elements 

Htri 

number  of  triangular  tubes  forming  the  Delta8  cell 

Nu 

Nusselt  number 

p 

total  or  partial  pressure  (N  m-2) 

Pel 

electrical  power  (DC)  (W) 

Qdot 

heat  flux  (Wm-1) 

r 

radius  (m) 

rj 

reaction  rate  of  reaction  j  [mols-1  m-1  ] 

rdl-reac 

diff.  limited  and  integration  length  specific  reaction 
rate  (mols-1  m-1) 

n 

inner  radius  (m) 

I'm 

middle  radius  (m) 

r0 

outer  radius  (m) 

ro 

reaction  order 

R 

resistance  (£2);  ideal  gas  constant  (J  mol-1 1<-1 ) 

T 

temperature  (K) 

fadtK 

gas  temperature  in  air  delivery  tube  (K) 

fal< 

anode  gas  temperature  (K) 

^amb 

ambient  temperature  (K) 

Tsk 

solid  structure  temperature  (I<) 

UF 

fuel  utilization 

Wc 

channel  width  of  planar  cell  channels  (m) 

Wce 

width  of  planar  cell  (m) 

Wq8 

Delta8  cell  width  (m) 

xp 

active  fraction  of  planar  cell  area  covered  by  bipolar 
rib 

xt 

inactive  fraction  of  tubular  cell  area 

y 

molar  fraction 

Greek  letters 

a 

upper  triangular  half  angle  (°) 

Wan 

convective  heat  exchange  coefficient  at  anode 

(Wm-2  K-1) 

Wca 

convective  heat  exchange  coefficient  at  cathode 
(W  m-2  I<-1 ) 

^insul 

heat  transfer  coefficient  through  insulation 
(W  m-2  I<-1 ) 

P 

transfer  coefficient 

A*, diff 

mass  transfer  coefficient  of  specie  x  (m  s-1 ) 

Ax.diff-reac  diff.  limited  reaction  conversion  coefficient  of 
specie  x 

Ax.reac 

reaction  conversion  coefficient  of  specie  x 

8 

thickness  of  component,  Delta8  cell  wall  thickness 

(m) 

AH 

heat  of  reaction  (J  mol-1 ) 

£ 

porosity  of  porous  media 

y 

pre-exponential  factor  for  exchange  current  density 

(Am”2) 

Pact 

activation  polarization  voltage  loss  (V) 

Pdiff 

diffusion  polarization  voltage  loss  (V) 

Pohm 

ohmic  voltage  loss  (V) 

X 

air-to-fuel  ratio 

hs 

solid  structure  heat  conductivity  coefficient 
(Wm-'K-1) 

Xan 

thermal  conductivity  of  anode  gas  (W  m-1  K-1 ) 

hca 

thermal  conductivity  of  cathode  gas  (W  m-1  K-1 ) 

Vav 

average  molecular  speed  (m  s-1 ) 

Vij 

stoiciometric  coefficient  of  specie  i  in  reaction  j 

P 

specific  conductivity  (1  m-1 ) 

X 

tortuosity  of  porous  media 

Subscripts  and  superscripts 

0,  in 

inlet 

act 

active 

ADT 

air  delivery  tube 

AEC 

anode-electrolyte-cathode  assembly 

an 

anode  gas  channel,  anode  electrode 

av 

average 

ca 

cathode  gas  channel,  cathode  electrode 

chact 

chemically  active 

circ 

circumferential 

conv 

convective 

cross 

cross-sectional 

cs 

catalyst  surface 

D8 

Delta8 

ed 

educts 

eff 

effective 

el 

electrolyte 

elact 

electrochemically  active 

equiv 

equivalent 

FC 

fuel  cell  tube 

hloss 

heat  loss 

hyd 

hydraulic 

ic 

interconnect 

ion 

ionization 

mix 

gas  mixture 

op 

operational 

oxi 

oxidation 

P 

planar 

prod 

products 

SH 

sensible  heat 

stoic 

stoiciometric 

STR 

steam  reforming 

t 

tubular 

tot 

total 

TPB 

triple  phase  boundary 

WGS 

water-gas-shift 
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models  in  systems  analysis  calculations  which  often  require  several 
hundred  iterative  steps  to  find  an  overall  solution  for  the  investi¬ 
gated  system.  Further,  CFD  models  require  an  exact  description  of 
the  cell  geometry  and  the  definition  of  the  calculation  mesh,  both 
being  time  consuming  tasks. 

We  propose  a  generalized  model  that  allows  for  the  fast  appli¬ 
cation  to  diverse  cell  designs,  considers  all  important  internal 
processes  and  at  the  same  time  features  short  calculation  times 
of  approximately  2  min  enabling  overall  system  simulations. 


2.  Model  definition 

Our  model  computes  temperature,  reactant  partial  pressure, 
current  density  and  voltage  loss  distributions.  In  addition,  integral 
values  such  as  the  absolute  power  output,  fuel  and  air  utilization 
are  predicted  for  a  user-defined  operational  voltage,  fuel  and  air 
input  flow  and  temperature.  Any  combination  of  H2,  CO,  C02,  H20, 
CFI4  and  N2  is  allowed  as  fuel  composition  and  the  oxygen  fraction 
in  the  cathode  gas  can  be  set  to  any  value.  The  electrochemical  con¬ 
version  of  hydrogen  and  carbon  monoxide  are  considered  to  occur 
in  parallel.  Methane  oxidation  through  reforming  reactions  is  also 
considered  via  a  Langmuir-FIinshelwood  kinetic  model. 

2.1.  Assumptions ,  definitions  and  model  structure 
The  main  assumptions  of  our  model  are: 

•  Steady-state  conditions  are  assumed. 

•  The  considered  SOFC  can  be  axially  divided  into  a  number  of  con¬ 
tinuously  stirred  tank  reactors  (CSTR),  termed  control  volumes 
(CV). 

•  Pressure  losses  in  the  gas  channels  are  not  considered. 

•  Due  to  their  high  electrical  conductivity,  the  interconnector  or 
bipolar  plates  of  planar  cell  designs  are  considered  equipotential 
[10]. 

•  The  flow  regime  in  the  gas  channels  is  assumed  as  plug  flow  [11  ]. 
Turbulent  flow  is  not  considered,  despite  its  possible  occurrence 
in  the  oxidant  gas  channels  due  to  excessive  flows  of  cooling  air 
[12].  Mass  transport  in  normal  direction  to  the  flow  direction  is 
assumed  to  be  by  diffusion  only. 

•  A  constant  Nusselt  number  with  the  value  4  is  assumed  [11  ]. 

•  Radiative  heat  transport  between  solid  structures  and  gas  streams 
is  ignored  [12]. 

•  The  heats  of  heterogeneous  reactions  are  attributed  to  the  fuel 
cell  solid  structure. 

•  The  diffusion  of  reactants  through  the  gas  channel  to  the  catalytic 
surface  is  considered. 

•  The  water-gas-shift  reaction  is  regarded  as  a  homogenous  reac¬ 
tion  at  equilibrium.  This  assumption  is  found  very  frequently  in 
the  literature,  e.g.  [10,13,14],  and  was  confirmed  in  experiments 
conducted  by  the  authors  (manuscript  in  preparation). 

Besides  local  values  such  as  solid  temperature  and  current  den¬ 
sity,  our  model  calculates  integral  values  which  characterize  the 
operational  conditions  of  the  investigated  cell  design  and  fuel 
composition.  The  definitions  are  given  in  the  following: 

•  The  fuel  utilization  is  defined  as  the  ratio  of  the  hydrogen  equiv¬ 
alent  mol  flows  at  the  cell  inlet  and  outlet. 


Our  model  consists  of  three  coupled  sub-models,  namely  the 
electrochemical  performance  model,  the  mass  balance  model  and 
the  energy  balance  model.  Fig.  1  depicts  the  model  structure  and 
the  iterative  solution  algorithm.  The  three  sub-models  are  coupled 
via  the  fuel  composition  and  the  temperatures  of  the  gases  and 
the  solid  structure  yielding  a  highly  non-linear  equation  system 
which  can  only  be  solved  iteratively.  In  this  work,  the  modeling 
package  “ATFIENA  Visual  Studio  11.0”  [15],  was  used  to  implement 
the  numerical  code  in  the  FORTRAN  language. 

The  electrochemical  performance  model  calculates  the  current 
density  for  a  given  fuel  composition  and  operational  voltage  in 
a  control  volume.  The  core  of  the  model  is  the  voltage  balance 
according  to  which  the  Nernst  voltage  minus  the  current  density- 
dependent  voltage  losses  has  to  equal  the  operational  voltage. 

The  mass  balance  model  requires  the  current  density  as  input 
for  the  calculation  of  the  conversion  rates  of  the  electrochemically 
active  species.  According  to  the  aforementioned  assumptions,  these 
are  hydrogen  and  carbon  monoxide.  Furthermore,  the  mass  balance 
model  includes  the  calculation  of  reaction  rates  for  all  considered 
homogenous  and  heterogeneous  reactions.  Besides,  the  mass  bal¬ 
ance  model  calculates  the  heat  source  terms  and  the  convective 
heat  transport  terms,  which  are  coupled  to  the  reactant  transport 
to  and  from  the  electrode  surfaces. 

The  energy  balance  model  serves  for  the  calculation  of  the  effec¬ 
tive  temperatures  of  the  solid  structure  and  the  gas  channels  based 
on  the  heat  source  terms  stemming  from  the  mass  balance  model. 
Furthermore,  the  energy  balance  model  includes  the  calculation  of 
the  purely  convective  heat  transfer  between  the  gases  and  the  solid 
structure  as  well  as  the  conductive  heat  transport  within  the  solid 
structure. 


2.2.  Gas  and  solid  material  properties 


2.2.1.  Thermodynamic  properties  of  gases 

The  calculation  of  the  thermodynamic  properties  of  the  involved 
species  is  the  basis  of  any  numerical  model  of  a  technical  process 
or  apparatus.  The  mean  temperature  and  the  temperature  varia¬ 
tion  occurring  in  the  process  determines  the  which  correlations 
are  the  most  suitable.  For  processes  which  are  operated  in  a  narrow 
temperature  band  it  is  valid  to  assume  temperature  independent 
material  properties.  This  approach  is  not  applicable  to  fuel  cell  mod¬ 
els  due  to  the  fuel  cell  operational  temperatures  which  can  vary 
from  700  °C  to  1000  °C,  see  [16]. 

For  the  estimation  of  the  heat  capacity,  thermal  conductivity 
and  viscosity  of  pure  gas  species,  correlations  according  to  the 
DIPPR  801  standard  (Design  Institute  for  Physical  Property  Data) 
were  implemented  as  subroutines  [17,18].  In  order  to  enable  the 
prediction  of  diffusion  limited  phenomena,  our  model  includes  a 
subroutine  for  the  calculation  of  effective  diffusion  coefficients. 
Two  cases  were  considered  to  describe  the  diffusion  of  gases 
through  the  electrodes  of  fuel  cells.  The  first  case  is  the  Knudsen 
diffusion,  which  is  of  importance  when  the  mean  free  path  length 
of  the  gas  molecules  is  large,  compared  to  the  pore  diameter.  In 
this  case,  the  molecules  collide  more  often  with  the  pore  walls  then 
with  each  other.  The  corresponding  Knudsen  diffusion  coefficient 
is  calculated  according  to: 


^av,* 


with  uav,x  = 


(2) 


UF  =  1  -  ?H2’°Ut  with  EqH 2  =  yHj  +  yco  +  4yCH4  ( 1 ) 

,  in 

The  hydrogen  equivalent  mol  flow,  Equ2,  is  calculated  assuming 
that  all  methane  is  completely  reformed  to  hydrogen  and  carbon 
dioxide. 


The  second  case  is  the  molecular  diffusion,  where  the  mean  free 
path  length  of  the  molecules  is  small  compared  to  the  pore  diame¬ 
ter.  In  this  case,  the  gas  molecules  collide  more  frequently  with  each 
other  than  with  the  pore  walls.  For  the  calculation  of  the  molecular 
diffusion  coefficients  the  method  according  to  the  kinetic  gas  the¬ 
ory  and  the  Fuller  method  were  implemented.  Detailed  information 
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Fig.  1.  Model  structure  and  iterative  solution  algorithm. 


about  both  methods  can  be  found  in  [19].  In  order  to  estimate  the 
molecular  diffusion  coefficient  of  a  specific  species  in  gas  mixtures, 
the  following  mixing  rule  was  implemented  [19]: 


D 


m,  mix,x 


l  -yx 

^  ^yj/Dm,i,x 


(3) 


2.2.2.  Temperature- dependent  solid  material  properties 

The  materials  used  in  solid  oxide  fuel  cells  have  to  fulfill  differ¬ 
ent  tasks  and  their  physical  properties  influence  the  behavior  of  the 
cell.  The  most  important  properties  in  this  respect  are  the  electrical 
and  the  ionic  conductivities,  which  have  a  direct  influence  on  the 
power  output  of  the  cell.  As  a  general  rule,  the  ceramic  materials 
employed  in  SOFCs  have  only  poor  electrical  and  ionic  conductivi¬ 
ties  at  room  temperatures.  At  elevated  temperatures  however  the 


conductivity  improves.  The  employed  correlations  for  the  temper¬ 
ature  dependent  specific  conductivity  were  proposed  in  [12]: 


Pan  — 


95.0  x  10b 


Tsk 


/— 1150.0\ 

exp(^H 


ln4  / 10300.0 
Pel  =  3-34  x  10  exp 


Pea  — 


42.0  x  10b 


Tsk 


■  exp 


V  TsK 

( 


-1200.0 


Pic  — 


9.3  x  10b 

?sK 


-  exp 


Tsk 


(4) 

(5) 

(6) 

(7) 
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2.3.  Geometries  and  discretization 

Fuel  cell  designs  under  development  mainly  differ  in  the  geome¬ 
try  and  the  employed  materials.  Concerning  the  geometry,  it  is  very 
useful  to  distinguish  between  the  macro-  and  the  microgeometry  of 
fuel  cells.  The  micro  geometry  describes  the  construction  of  the  AEC 
assembly,  while  the  macrogeometry  describes  the  construction  of 
the  cell  itself. 

The  micro  geometry  has  a  strong  influence  on  the  electrochem¬ 
ical  performance  of  the  fuel  cell  via  ohmic  losses  and  transport 
limitations  due  to  the  diffusion  of  the  reacting  species  through  the 
electrodes  to  the  sites  where  the  electrochemical  reactions  occur. 

The  macrogeometry  has  an  impact  on  the  electrochemical  per¬ 
formance  via  ohmic  losses  occurring  in  the  current  conducting 
components  of  the  fuel  cell,  the  heat  balance  via  convective  and 
conductive  heat  exchange  areas  and  the  mass  balance  via  the  cat- 
alytically  active  surfaces  areas. 

Accounting  for  the  strong  influence  of  the  micro-  and  macro¬ 
geometries  of  different  fuel  cell  designs  on  the  cell  characteristics 
was  the  most  important  requirement  for  our  model.  The  chosen 
approach  for  describing  the  geometry  is  based  on  characteristic 
lengths  and  areas.  The  corresponding  definitions  for  the  considered 
geometries  are  presented  in  the  following. 


structures  according  to  the  number  of  gas  channels.  The  follow¬ 
ing  characteristic  lengths  and  areas  were  defined  for  the  repeating 
structure  of  the  planar  geometry: 


•  Total  active  area,  AP(act 
ap’act=® ,c 

•  Electrochemically  active  length,  Jp  elact 

fp.elact  =  Wc  +  (  (^)  -  Wc)  Xp.elact 

•  Chemically  active  length,  /p>Chact 

fp,chact  —  wc  +  (  (  n~h  )  _  WC)  Xp’chact 

•  Hydraulic  diameter  of  anode  gas  channel,  dp,hyd,an 

A  _  2wchPian 

ap,hyd,an  -  ,  u 

Wc  +  Up, an 

•  Hydraulic  diameter  of  cathode  gas  channel,  dPihyd,ca 


(8) 


(9) 


(10) 


(ID 


2.3 A.  Planar  geometry 

Fig.  2  depicts  the  considered  planar  geometry.  The  planar  cell 
stack  consists  of  bipolar  plates  which  form  the  rectangular  gas 
channels.  The  bipolar  plates  are  separated  by  AEC  assemblies.  The 
surface  of  the  AEC  assemblies  is  consequently  partially  covered  by 
the  ribs  of  the  bipolar  plates.  The  planar  cell  is  divided  into  repeating 


d 


p,hyd,ca 


2wcfrp,ca 

Wc  +  flp.ca 


(12) 


•  Circumferential  length  of  anode  channel  perpendicular  to  gas 
flow,  (p,circ,an 


Jp,circ,an  —  2(wc  +  hp.an) 


(13) 


Interconnect^/ 
£5^ 


Cathode  channel 


Anode  channel  -T 
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Cathode  channel^^j" 


Anode  channel-p^ 


Macro  and  micro 
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repeating  structure 
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on  anode  side 


j  Thickness  of  anode 
Thickness  of  electrolyte 
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Interconnect 
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cathode  side 


Fig.  2.  Generalized  planar  cell  geometry  and  repeating  structure  description. 
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•  Circumferential  length  of  cathode  channel  perpendicular  to  gas 
flow,  (PiCirCiCa 

fp,circ,ca  =  2(wc  +  Jlp.ca)  (14) 

•  Cross-sectional  area  of  solid  structure  perpendicular  to  gas  flow 
direction,  Across 


^p,  cross 


(WCe  \ 

J  (^ic,an  +  ^p,an  +  <5p,el  +  ^p,ca  +  ^ic,ca) 


— Wc(hp?an  +  frp,ca) 


(15) 


2.3.2.  Tubular  geometry 

Fig.  3  depicts  the  standard  tubular  geometry.  The  standard  tubu¬ 
lar  cell  consists  of  a  AEC  tube  with  one  closed  end  and  a  concentrical 
air  delivery  tube  (ADT).  The  annular  gap  between  the  outer  surface 
of  the  ADT  and  inner  surface  of  the  AEC  tube  forms  the  cathode  gas 
channel.  The  anode  gas  flows  around  the  outer  surface  of  the  tube. 
Thus,  the  geometry  of  the  anode  gas  channel  is  not  defined  by  the 
cell  design  itself  but  rather  by  the  arrangement  of  the  single  tubes  in 
the  stack.  For  modeling  purposes  however,  an  annular  gap  shaped 
anode  gas  channel  was  assumed.  The  fuel  cell  tube  itself  was  consid¬ 
ered  as  a  repeating  structure.  The  following  characteristic  lengths 
and  areas  were  defined  for  the  tubular  geometry: 

•  Middle  radius  of  the  AEC  assembly,  rmAEC 

^  ri,FC  +  (^t,ca  +  ^t,el  +  ^t,an)  ,  . 

rm,  AFC  =  - 2 -  (lb) 

•  Outer  radius  of  the  cell  tube,  r0,Fc 

ro,FC  —  ri,FC  +  ^t,ca  +  $t,el  +  ^t,an  (17) 

•  Total  active  area,  At,act 


At, act  =  27rrm?AEC/t  (18) 

•  Electrochemically  active  circumferential  length,  /t  elact 

^t,elact  —  27rrm?AEC(l  —  *t,ic)  (19) 

•  Chemically  active  circumferential  length,  ltchact 

^t,chact  —  27Tr0?Fc(l  ^t,  ic )  (20) 

•  Elydraulic  diameter  of  air  delivery  tube,  dt,hyd,ADT 

dt,hyd,ADT  =  2ri?ADT  (21) 


•  Elydraulic  diameter  of  virtual  anode  gas  channel,  dt,hyd,an 

^t,hyd,an  =  2  •  fat,  an  (22) 

•  Elydraulic  diameter  of  cathode  gas  channel,  dt,hyd,ca 

dt,hyd,ca  =  2(ri  FC  -  r0?ADT)  (23) 

•  Circumferential  length  of  air  delivery  tube  perpendicular  to  gas 
flow,  /tiCirc,ADT 

^t,circ,ADT  =  7t(ro,ADT  +  ri,ADl)  (24) 

•  Circumferential  length  of  virtual  anode  channel  perpendicular  to 
gas  flow,  /t)CirCian 

^t,circ,an  —  27Tr0,Fc(l  —  *t,ic)  (25) 

•  Circumferential  length  of  cathode  channel  perpendicular  to  gas 
flow,  /t(CirCiCa 

^t,circ,ca  —  ^Ttr^  Ec  (26) 

•  Cross-sectional  area  of  solid  structure  perpendicular  to  gas  flow 
direction,  Across 

cross  =  (7T rji;  fc)  -  (7Tr?FC)  (27) 

2.3.3.  Triangular  tube  geometry 

Fig.  4  depicts  the  triangular  tube  geometry  (D8).  It  is  assumed 
that  similar  to  the  standard  tubular  cells,  the  D8  employs  centered 
air  delivery  tubes  (ADT)  for  the  injection  of  the  cathode  air.  The 
D8  cell  is  made  up  of  eight  connected  triangular-shaped  AEC  tubes 
with  triangular  cross-section  and  closed  ends,  Fig.  4a.  The  spaces 
between  the  ADTs  and  inner  surfaces  of  the  triangular  AEC  tubes 
form  the  cathode  gas  channels,  Fig.  4d.  The  outer  upper  surface  of 
the  D8  cell  and  the  outer  lower  surface  of  the  next  D8  cell  form  the 
anode  gas  channels,  Fig.  4d.  Note  that  the  D8  cells  are  stacked  like 
planar  cells,  Fig.  4a.  For  better  contacting  and  current  distribution, 
a  porous  contact  layer  is  placed  between  the  individual  cells.  Due 
to  its  high  electrical  conductivity  and  porosity,  the  contact  layer 
was  not  considered  in  our  model.  The  D8  repeating  structure  was 
defined  as  one  triangular-shaped  AEC  tube. 

The  description  of  the  triangular  cell  design  macro  geometry 
is  more  complex  than  that  of  the  standard  tubular  and  planar  cell 
design.  However,  the  complexity  can  be  reduced  by  introducing 
geometrical  dependencies  for  the  curvature  radii.  For  our  model, 
it  was  assumed  that  the  upper  outer  curvature  radius  r0  equals 
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Fig.  4.  Triangular  tube  cell  geometry. 


1.5  times  the  D8  cell  tube  wall  thickness  8.  The  lower  outer  cur¬ 
vature  radius  rm  was  assumed  equal  to  the  cell  tube  wall  thickness 
and  the  inner  curvature  radius  rx  equal  to  half  the  cell  tube  wall 
thickness,  Fig.  4b.  With  the  mentioned  assumptions  the  follow¬ 
ing  characteristic  lengths  and  areas  were  derived  for  the  triangular 
geometry: 

•  Electrochemically  and  chemically  active  width,  /ps.eiact 
28  90  -  a 


D8,  elact  —  ^D8,  chact  — 


+  ■ 


x  = 


ntricosa  ntri90 
wDg  5  8  90  -  oi 


7x8 +  2X  with 


ntr i  2  tana 


+  ■ 


72 


-n8 


(28) 


The  active  width  of  one  repeat  element  is  the  arithmetic  average 
of  the  overall  active  width  of  the  D8  cell.  Note  that  the  additional 
electrochemically  active  areas  of  the  lower  curvature  of  the  two 
outer  triangular  cells  were  hence  considered,  Fig.  4d. 

Total  active  area,  AD8>act 


^D8,act  —  ^D8,elact^D8 


(29) 


Circumferential  length  of  ADT  perpendicular  to  gas  flow, 

^D8,circ,ADT 


^D8,circ,ADT  =  ^(*08,0,  ADT  +  ^D8,i,ADT ) 


(30) 


Circumferential  length  of  anode  channel  perpendicular  to  gas 
flow,  /D8,circ,an 


(31) 


Circumferential  length  of  cathode  channel  perpendicular  to  gas 
flow,  /D8,circ,ca 


D8,circ,ca 


sin  a  +  1  (  wD8  8  +  2  sin  a 


sin  a 
90  +  a 


Htri 


cos  a 


7x8  - 


38 


120  tan  (45  -  f ) 

1  Cross-sectional  area  of  cathode  channel,  AD8i, 
.  X2 

Ad8,  cross. 


(32) 


4  tan  a 


3  „2  /  90  +  a 

+  (  — -t^tx  -  cot 


360 


-nr. 


D8,o,ADT 


with  X= 


(«-?)) 

wD8  5 +  2  sin  a 


^tri 


cos  a 


(33) 


Cross-sectional  area  of  anode  channel,  AD8>cr0ss,an 
a  3  o7  90  —  a  82 

A™’ cross, an  =  X2tana  +  yiS2  -  ^  -  — 

x  =  WD8  cos«  +  nt|-i3(sin«-2)  and 


with 


Y  = 


2 ntri  sin  a 
2sina  +  cosa  -  1  -  sin2  a 


sin  a  cos  a 

Cross-sectional  area  of  solid  structure,  AD8  c 
WD8 


A)8,c 


x  = 


nt 


-X  —  A[)8  cross  ca  —  Ap8 


cross,  an 


-  7rr; 


D8,o,ADT 


(34) 


with 


WD8  cos  a _ 8_  3  ^ 

2ntrisina  sin  a  2 


/  _  /  ,  Wp8 

*D8,circ,an  —  *D8,elact  ' 

“tri 


(35) 
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Fig.  5.  Equivalent  circuit  of  the  electrochemical  performance  model. 


Hydraulic  diameter  of  ADT,  dD8,hyd,ADT 

^D8,  hyd,  ADT  =  ^rD8,  i,  ADT 

Hydraulic  diameter  of  anode  gas  channel,  dD8hyda 

4AD8  ,  cross,  an 


uD8,hyd,an  : 


D8,cric,an 


Hydraulic  diameter  of  cathode  gas  channel,  dD8hydca 

4Ad8,  cross,  ca 


uD8,hyd,ca  — 


^D8,cric,an  +  2jrrD8,i,  ADT 


(36) 


(37) 


(38) 


Hence,  the  model  must  distinguish  currents  originating  from  the 
oxidation  of  hydrogen,  /h2,  or  of  carbon  monoxide,  /c o,  which  sum 
up  to  the  total  produced  current,  hot- 


hot  =  I  h2  +  ho 


(42) 


From  Eq.  (41 )  and  Fig.  5  it  can  be  seen,  that  a  current  originating 
from,  e.g.  carbon  monoxide  oxidation  can  only  be  produced,  when 
the  according  Nernst  voltage  minus  the  voltage  losses  of  the  oxygen 
electrochemistry  and  the  ohmic  losses  due  to  an  already  flowing 
current  originating  from  hydrogen  oxidation  stills  yields  a  positive 
value. 


2.4.  Electrochemical  performance  model 

The  electrochemical  performance  model  calculates  the  current 
density  for  a  given  fuel  composition  and  operational  voltage  con¬ 
sidering  activation,  ohmic  and  diffusion  losses.  Fig.  5  shows  the 
equivalent  circuit  of  the  electrochemical  performance  model. 

The  starting  point  for  the  calculation  of  the  current  density  is  the 
reversible  potential  of  the  fuel  gas,  referred  to  as  the  Nernst  voltage. 
The  Nernst  voltage  of  the  hydrogen  oxidation,  Nernst,  h2  *  and  carbon 
monoxide  oxidation,  Nernst, co>  is  calculated  using  Eq.  (39),  where 
n  stands  for  the  number  of  exchanged  electrons  per  reaction.  For 
the  oxidation  of  hydrogen  or  carbon  monoxide  it  Eq.  (2),  while  it 
equals  4  for  the  ionization  of  oxygen. 

rTsk  \ 
nF  ) 

i  =  H2,  CO  and  j  =  H20,  C02  (39) 

The  values  of  the  equilibrium  constants  for  the  hydrogen  and  carbon 
monoxide  oxidation  reaction  were  taken  from  [20].  The  tabulated 
values  were  compiled  to  yield  two  fit  correlations,  which  both  have 
the  form  of  Eq.  (40): 

X=y0+^i  exp  (-^)  +A2 exp  (-^)  +A3exp  (-^)  (40) 

X  denotes  either  the  equilibrium  constant  of  the  hydrogen  oxidation 
or  of  the  carbon  monoxide  oxidation,  depending  on  the  employed 
set  of  fit  correlation  coefficients,  Table  1.  Fig.  5  depicts,  that  the 
Nernst  voltages  of  the  hydrogen  and  carbon  monoxide  oxidation 
minus  the  current  density-dependent  voltage  losses  have  to  equal 
the  operational  voltage  of  the  fuel  cell  (Eop): 

£0p  =  ^Nernst,  i  —  4act,i  —  4diff,i  —  4act,02  —  4diff,02  —  4ohm  with 

i  =  H2,  CO  (41) 


^Nernstq  : 


ln/Cpioxi  +  In 


PiP ; 


02 


Pi 


with 


2.4.1.  Activation  polarization 

Even  when  no  net  current  is  drawn  from  the  fuel  cell,  the 
electrochemical  reactions,  Eqs.  (43)  and  (44)  (anode)  and  Eq.  (45) 
(cathode),  occur  at  equal  rates  in  both  directions: 


H2+02-+*  H20  +  2e~ 
CO  +  02“  C02  +  2e~ 
02  +  4e~  4*  202~ 


(43) 

(44) 

(45) 


The  so-called  exchange  current  represents  the  current  flowing  in 
both  directions  at  equilibrium  conditions.  In  order  to  generate  a 
net  current  into  either  direction  a  certain  potential  needs  to  be 
applied.  This  potential  is  referred  to  as  activation  polarization,  pact, 
and  mathematically  best  described  by  the  Butler-Volmer  equation 
[21]: 


I  =  Io 


exp 


(onFr/za 


) 


-  exp 


(-P-/8) 


nFr)  act 
RTsK 


) 


(46) 


The  Butler-Volmer  equation  comprises  two  important  model 
parameters.  The  first  parameter  is  the  transfer  coefficient  which 
represents  the  part  of  the  change  in  polarization  leading  to  a  change 
in  the  reaction  rate  constant.  For  fuel  cell  applications  and  in  our 


Table  1 

Coefficients  of  the  equilibrium  constant  fit  correlation  for  hydrogen  and  carbon 
monoxide  oxidation 
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Table  2 

Equivalent  resistance  formulas  for  five  cell  sub-units  with  the  arrowed  lines  indicating  the  current  flow  path  from  cathode  to  anode  [25] 
Cell  sub-unit 


Geometry  factor  J 


Resistance  R(Qm) 


s 


if 

if 

p 

P 

p 

X 

ft 

X 

,\v 

l\X 

pi 

8 


2^5 


Pc  a/ 8  c  a 
Pel^el 


Pel^el(Pca/8  ca) 

tanhj 


I  1 


ESS 


::  i 


5X3 


V?i(£+e) 


((Pan/ ^an)2  +  (Pea/ <5ca)2)  ■  COShJ  +  (panPca/ <5an<5ca)  ■  (2  +J  ■  sinhj) 
'sj ( VPel^el) '  \J (Pan/ ^an)  +  (Pea/ <5ca)  ■  sinhj 


Pan  <5  an  +  Pel^el  +  Pea  <5  ca 


L\l — 

V  Pel  "el  \  ^an  ^ca  / 


Pel^el((Pan/ <$an)  +  (Pea/ <5ca)) 


tanhj 


model,  the  value  is  assumed  to  be  0.5  [21,22].  The  second  parameter 
is  the  exchange  current  density  /0,  which  represents  the  forward 
and  reverse  electrode  reaction  rate  at  the  equilibrium  potential. 
High  exchange  current  densities  mean  a  high  electrochemical  reac¬ 
tion  rate  and  low  activation  losses  [23].  A  large  number  of  models 
to  calculate  the  exchange  current  density  at  the  anode  and  cath¬ 
ode  can  be  found  in  the  literature.  In  the  present  work,  the  anode 
exchange  current  density  is  calculated  using  Eq.  (47)  for  the  hydro¬ 
gen  and  the  carbon  monoxide  electrochemistry.  The  exponents  of 
the  molar  fractions  near  to  the  stoichiometry  of  the  anodic  elec¬ 
trochemical  reactions  result  in  a  direct  proportional  dependence  of 
the  anode  exchange  current  density  to  the  educt  and  an  inversely 
proportional  dependence  to  the  product  partial  pressures.  With 
the  inverse  dependency  of  the  activation  losses  to  the  exchange 
current  density,  the  exponents  of  the  molar  fractions  yield  lower 
activation  losses  with  increased  educt  or  decreased  product  par¬ 
tial  pressures  and  vice  versa  [22].  The  cathode  exchange  current 
density  is  given  by  Eq.  (48).  The  exponent  of  the  molar  fraction 
yields  lower  cathode  activation  losses  with  increased  oxygen  partial 
pressures. 

io,,i  wi,h 

i  =  H2,  CO  and  j  =  H20,  C02  (47) 


It  has  to  be  emphasized,  that  due  to  the  implicit  character  of 
the  Butler-Volmer  equation  its  solution  can  only  be  determined 
numerically.  This  effort  is  however  justified  by  higher  accuracy 
compared  to  simplified  equations  such  as  the  Tafel  equation  [24]. 

2.4.2.  Ohmic  polarization 

Electronic  currents  through  the  electrodes  and  interconnects  as 
well  as  ionic  currents  through  the  electrolyte  induce  frictional  heat 
and  thereby  voltage  losses.  The  calculation  of  the  ohmic  voltage 
losses  simply  follows  Ohm’s  law: 

4ohm  —  ^equiv^tot  (49) 

The  equivalent  resistance  (Re quiv)  of  SOFCs  depends  on  the  geom¬ 
etry  and  the  conductivity  of  the  current  conducting  components; 
see  Eqs.  (4)-(7).  In  the  present  model,  the  area  specific  equivalent 
resistance  (ASR)  is  approximated  using  the  transmission  line  model 
instead  of  calculating  it  via  the  numerical  solution  of  the  Laplace 
equation.  The  methodology  was  developed  by  Nisancioglu  [25]  and 
others  [26],  and  features  a  high  flexibility  towards  different  fuel 
cell  designs.  The  general  idea  of  the  transmission  line  model  is  that 
complex  geometries  can  be  described  by  a  set  of  five  standard  cell 
sub-units.  Each  of  these  cell  sub-units  has  a  characteristic  current 
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path  pattern  and  according  analytical  formula  for  the  equivalent 
resistance  calculation,  Table  2. 

2.4.2.L  Application  to  tubular  geometry.  Fig.  6  shows  the  cross  sec¬ 
tion  with  current  paths  of  the  tubular  geometry  considered  in  our 
model,  Fig.  6a.  The  transmission  line  model  approximation,  Fig.  6b, 
consists  of  two  type  3  cell  sub-units  connected  in  parallel,  repre¬ 
senting  the  two  halves  of  the  AEC  tube,  and  one  type  2  cell  sub-unit, 
representing  the  interconnect  area,  connected  in  series.  Due  to  the 
perfect  symmetry  of  the  tubular  geometry,  the  equivalent  resis¬ 
tance  is  only  calculated  for  one  half-cell  and  then  divided  by  two. 
The  length  of  the  type  3  cell  sub-unit  (Le)  for  one  half-cell  is  cal¬ 
culated  according  to  Eq.  (50)  and  for  the  cell  sub-unit  type  2  (Lic) 


according  to  Eq.  (51): 

Le  =  ^m,AEc(l  —  jc) 

(50) 

Lie  =  ^t^m,AEC^t,ic 

(51) 

The  resistance  of  one  half  of  the  cell  tube  (Re)  is  then  given  by 

((Pan/^an)  +  (Pca/<^ca)  )  COShJe 
n  +(PanPca/^an^ca)(2 +Je  Sinhje)  .  u 

Rc  =  — with 
\j ( 1/Pel^el)  \J (Pan/ ^an)  +  (Pea/ $ca)  SinhJe 


The  resistance  of  the  interconnect  ( Ric )  is  calculated  according  to 

Ric  =  With  Jic  =  Lic Va/PdSdXPc/Sc)  ( 53 ) 

The  equivalent  resistance  is  then  calculated  following  Kirchhoff  s 
law  for  in-series  connection.  In  order  to  obtain  the  ASR  value  for 
the  calculation  of  the  ohmic  voltage  loss  using  a  current  density, 
the  equivalent  resistance  is  multiplied  with  the  electrochemically 
active  circumferential  length  of  the  cell  tube: 

^t,equiv  —  ^t,elact(^e  +  ^ic)  (54) 


2A.2.2.  Application  to  triangular  geometry.  The  triangular  geome¬ 
try  exhibits  current  paths  which  can  be  described  using  the  same 
approach  as  for  the  standard  tubular  cell,  Fig.  7.  The  cell  sub-unit 
lengths,  Le  and  Lic,  are  given  by 


^D8,e  —  0-5/D8,elact 


(55) 


^D8,ic 


Wps 

2 n^i 


(56) 


The  ASR  is  obtained  by  multiplying  the  equivalent  resistance 
with  half  of  the  electrochemically  active  length  of  the  repeat 
element: 


^D8,equiv  —  ^.5/pg  eiact(Re  +^ic) 


(57) 


Fig.  7.  Cross  section  of  the  triangular  tube  cell  design  (a)  and  transmission  line  model  (b). 
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2.423.  Application  to  planar  geometry.  Fig.  8  shows  a  cut  out  of 
the  considered  planar  stack  and  the  cross-section  of  a  repeating 
structure  with  current  paths,  Fig.  8a,  and  transmission  line  model 
approximation,  Fig.  8b.  The  equivalent  resistance  of  the  repeat  ele¬ 
ment  includes  six  cell  sub-units.  In  our  model,  the  current  flow 
through  the  AEC  assembly  is  divided  into  a  sub-unit  type  4  ( R i )  and 
type  5  (R2): 

p  2(pCa^ca  +  Pel^el  +  Pan^an)  rc-o\ 

wc/(nch-wch) 


R2 


\J  Pel^eMP  an/ $an  )  +  (Pca/<U) 

tanhj2 


with 


The  two  cell  sub-units  are  connected  in  parallel.  The  interconnect 
contacts  anode  and  cathode  of  two  adjacent  AEC  assemblies.  The 
resistances  of  the  interconnect  ribs  on  the  anode  (R3)  and  cathode 
(ft6)  side  are  calculated  using  type  1  cell  sub-units: 


Rs 


Re 


2pichp?an 
(wc/nch)  -  wch 

2pic/lp,ca 
(wc/nch)  -  wch 


(60) 

(61) 


The  resistance  of  the  bar  between  adjoining  ribs  is  calculated  by 
assuming  a  parallel  connection  of  a  type  1  cell  sub-unit  (R4)  and  a 
modified  type  5*  cell  sub-unit  (R5)  [27]: 


R4 


2pic((hic,an  +  fojc.ca)  ~  (^p,an  +  hp,ca)) 
(wc/nch)-wch 


(62) 


2.4.3.  Diffusion  polarization 

The  diffusion  polarization  results  from  the  difference  of  the 
Nernst  voltage  calculated  for  the  reactant  partial  pressures  in  the 
bulk  gas  phase  and  at  the  reaction  sites  located  at  the  triple  phase 
boundary  (TPB)  where  electrolyte  and  electrode  material  as  well  as 
reactant  gases  meet.  The  partial  pressures  of  the  educts  are  gener¬ 
ally  lower  while  those  of  the  products  are  generally  higher  at  the 
TPB  than  in  the  bulk  gas  phase.  The  exact  partial  pressure  values 
depend  on  the  local  current  density  and  the  material  parameters  of 
the  porous  electrodes  through  which  the  educts  and  products  have 
to  diffuse  in  order  to  get  to  the  TPB.  The  current  density  determines 
the  concentration  gradient  between  the  TPB  and  the  bulk  gas  phase. 
The  material  parameters  of  the  porous  electrodes  include  tortuos¬ 
ity,  porosity  and  pore  radius.  These  material  parameters  together 
with  the  fuel  composition  determine  the  effective  diffusion  coef¬ 
ficient  which  is  required  to  calculate  the  partial  pressures  at  the 
TPB  using  e.g.  Fields  law  of  diffusion.  Concentration  polarization 
becomes  important  for  highly  diluted  fuel  gases  and  at  high  current 
densities,  where  they  increase  against  an  asymptotic  maximum.  At 
this  specific  point  a  further  increase  of  the  current  is  impossible 
since  the  reactant  partial  pressures  at  the  TPB  equal  zero  due  to 
instant  conversion. 

Our  model  considers  diffusion  limitation  in  the  anode  and  the 
cathode.  At  the  anode,  the  educts  and  the  products  of  the  electro¬ 
chemical  reactions  diffuse  with  equal  rates  to  the  TPB  and  back 
into  the  gas  flow  channel  respectively.  This  equimolar  counter¬ 
flow  diffusion  yields  a  zero  net  diffusion  flux.  Applying  Fields 
law  of  diffusion  by  assuming,  that  the  electrochemical  reaction 
rate  equals  the  mass  flux  via  diffusion,  the  partial  pressures  of 
educts  and  products  at  the  TPB  can  be  straightforwardly  calcu¬ 
lated  based  on  the  partial  pressures  in  the  bulk  gas  phase  according 
to 


Rs 


_ Pic _ 

0.41(1  —  exp(  — 0.6(wch/((hican  +  h[c<ca)  —  (hp,an  +  hp,ca))))) 

(63) 


Finally,  the  equivalent  circuit  depicted  in  Fig.  8  yields  Eq. 
(64).  Similar  to  the  other  geometries,  the  ASR  is  multi¬ 
plied  with  the  electrochemically  active  length  of  the  repeat 
structure. 


R 


p,eqiv_ 


1  1 


y'm 


1  1 

+ 


(64) 


P™  =  Pi  T 


f  RTaK^an^an 
y  nFDe ff£an 


with  -for  i  =  H2,  CO  and 


+  fori  =  FI20,  C02 


(65) 


The  effective  diffusion  coefficients,  Deff)H  and  Deff  co,  are  calculated 
based  on  binary  molecular  diffusion  coefficients  of  hydrogen  and 
carbon  monoxide  with  all  other  considered  species  in  the  gas  mix¬ 
ture  and  the  Knudsen  diffusion  coefficient.  The  binary  molecular 
diffusion  coefficients  are  summed  up  to  yield  the  effective  molec¬ 
ular  diffusion  coefficient  according  to  Eq.  (3).  For  the  calculation 


154 


F.P.  Nagel  et  al.  /  Journal  of  Power  Sources  184  (2008)  143-164 


of  the  effective  diffusion  coefficient  it  is  assumed,  that  Knudsen 
and  molecular  diffusion  occur  in  parallel.  Therefore,  the  Bosanquet 
formula  was  implemented  in  the  present  model: 


^eff.i  — 


^m,i-mix  '  ^K,i 
^m,i-mix  +  ^I<,  i 


with  i  =  H2 ,  CO 


(66) 


Knowing  the  partial  pressures  of  the  educts  and  products  in  the 
bulk  gas  phase  and  at  the  TPB,  the  diffusion  voltage  losses  can  be 
computed  according  to 


with  i  =  H2,  CO  and 


j  =  H20,  C02 


(67) 


Table  3 

Coefficients  of  the  equilibrium  constant  fit  correlation  for  water-gas-shift  and 
methane  steam  reforming  reactions 


The  reaction  rates  of  the  WGS  and  REV-WGS  are  computed  via  Eqs. 
(75)  and  (76): 


At  the  cathode,  only  oxygen  takes  part  in  the  electrochemical  reac¬ 
tions,  where  it  is  consumed  without  forming  products. 


PToP2B-P 


A<,0 2(aP™  -  P)  -  An,02-N2 
C>K,02(^P02  ~P)  -  An,02-N2 


=  exp(/totTci<B)  with 


A  = 


B  = 


Po2  -P 


and 


TCaP£ca(£>m,02-N2  -  A<,02(^  -  1))  \  ^ 0 2 


eca2nFDm>o2_N2DK>o2  p 


Mn 


(68) 


Eq.  (68)  was  implemented  in  our  model  for  the  implicit  calcula¬ 
tion  of  the  partial  pressure  of  oxygen  at  the  cathode  TPB  [28].  In 
contrast  to  the  diffusion  at  the  anode,  the  effective  diffusion  coeffi¬ 
cient  is  calculated  for  the  non-dilute  two-component  gas  mixture  at 
the  cathode  according  to  [29]  which  includes  the  relation  between 
the  molecular  masses,  a,  of  the  components  nitrogen  and  oxygen. 
The  diffusion  voltage  loss  at  the  cathode  electrode  is  calculated 
analogously  to  the  anode  electrode: 


rWGs  =  10000.0  mol  m  1  s  1  bara  2  PcoPh2o  (75) 

tREV-WGS  =  10000.0  mol  m-1  s-1  bara-2  Pco2Ph2 

PcoPh2o  \  ,76n 

J<PwgsPco2Ph2  ) 

Table  3  lists  the  coefficients  for  the  equilibrium  constant  fit  corre¬ 
lation,  Eq.  (40),  of  the  WGS.  The  equilibrium  constant  values  were 
taken  from  [20].  In  contrast  to  the  WGS,  the  STR  was  considered 
as  heterogeneous  reaction.  An  applied  kinetic  model,  Eq.  (77),  was 
implemented  into  the  model.  The  implemented  applied  kinetic 
model  for  the  methane  steam  reforming  was  developed  by  Drescher 
[30],  who  related  the  reaction  rate  to  the  mass  of  catalyst.  Conse¬ 
quently  the  pre-exponential  factor  had  to  be  recalculated  to  an  area 
specific  value  in  order  to  allow  for  the  application  in  the  present 
model.  We  therefore  assumed  that  only  the  upper  25  microns  of 
the  nickel-cermet  participate  in  the  steam  reforming  reactions.  All 
other  data  for  the  recalculation  were  taken  from  [30]  in  order  to 
keep  the  applied  kinetic  model  concise: 


Pdiff,02  = 


2.5.  Mass  balance  model 


(69) 


rCH4STR,Dre 


288.52  mol  m-2  s-1  bara-2  pCH4  Ph2o 
xexp(-11000J  mol-1 /RTsK) 

1  +  16.0 1  bara-1  Pch4  +  0.143 1  bara-1  pn2o 
x  exp(39000J  mol-1  /RTsK) 


The  current  density  values  computed  by  the  electrochemical 
performance  model  are  related  to  the  reactions  Eqs.  (43)-(45)  via 
the  Faraday  law  yielding  the  area  specific  reaction  rates  of  the  elec¬ 
trochemical  reactions.  In  the  present  model,  the  axial  length  of  the 
gas  channels  was  defined  as  the  integration  variable.  In  order  to 
account  for  the  two-dimensional  distribution  of  the  electrochem¬ 
ical  reactions,  the  area  specific  reaction  rate  has  to  be  multiplied 
with  the  electrochemically  active  length  of  the  considered  fuel  cell 
design,  Zelact,  so  as  to  be  converted  to  the  integration  length  specific 
reaction  rate.  Accordingly,  Eqs.  (70)  and  (71)  were  applied  in  this 


model: 

ri,oxi  —  ^eiact  with  i  =  H2,  CO  (70) 

ro2ion  =  Zelact  (^)  (71) 

Besides  the  mentioned  electrochemical  reactions,  this  model  con¬ 
siders  water-gas-shift,  Eqs.  (72)  and  (73),  and  methane  steam 
reforming  reactions: 

CO  +  H20  C02  +  H2  (72) 

C02  +  H2*>  CO  +  H20  (73) 

CH4  +  H20  ^  CO  +  3H2  (74) 


Depending  on  the  geometry  of  the  anode  gas  channel,  the  diffusion 
of  the  reacting  species  from  the  bulk  gas  phase  to  the  catalyst  sur¬ 
face  can  be  slower  than  the  actual  chemical  reaction.  Therefore,  the 
partial  pressures  of  the  reactants  at  the  catalyst  surface  can  differ 
significantly  from  the  bulk  gas  phase.  The  calculation  of  the  exact 
reactant  partial  pressures  at  the  catalyst  surface  requires  the  solu¬ 
tion  of  the  complete  concentration  field  perpendicular  to  the  gas 
flow.  In  order  to  avoid  solving  the  corresponding  partial  differen¬ 
tial  equations,  a  mass  transfer  analogy  proposed  in  [11]  was  used 
to  compute  the  reactant  partial  pressure  at  the  catalyst  surface: 


_  Aqdiff  /  _  c 

Gqdiff  ^  \Px  Px 


(78) 


The  mass  transfer  coefficient  )diff  is  calculated  by  considering  the 
analogy  of  heat  and  mass  transfer  according  to 


A*,diff  — 


Pm, mix, x  m 

“t/p,hyd,an 


(79) 


The  calculation  of  the  reactant  partial  pressure  at  the  catalyst  sur¬ 
face  can  be  avoided,  by  considering  the  diffusion  process  and  the 
chemical  reaction  as  two  processes  connected  in  parallel.  Mathe¬ 
matically  this  requires  the  introduction  of  a  conversion  coefficient 
for  the  considered  chemical  reaction,  px,reac  which  is  based  on 
the  reaction  rate,  rYeac,  calculated  with  the  bulk  gas  phase  partial 
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pressures,  Eq.  (80).  Note  that  the  reaction  rate  is  divided  by  the  par¬ 
tial  pressure  of  the  diffusion-limited  specie,  px,  to  the  power  of  its 
reaction  order,  rox: 


6  -  r  RTsK 

Px,  reac  -  'reac  — fo^ 
Vx 


(80) 


The  diffusion  limited  and  integration  length  specific  reaction  rate, 
tdi-reac*  is  obtained  from 


Px0x 

hdl-reac  =  ^chactA>c,diff-rea  cpjr~^ 


With  ^,diff_reac  = 


Ae,diffA*,reac 
fix,  diff  +  Px,  reac 

(81) 


The  partial  pressure  of  the  diffusion-limited  species,  px,  to  the 
power  of  its  reaction  order,  rox,  is  thus  reintroduced.  Knowing  the 
reaction  rate  of  all  reactions  and  neglecting  axial  diffusion  mass 
transport,  the  spatial  distribution  of  the  species  along  the  anode 
channel  can  be  computed  according  to  the  differential  equation: 


^  withi  =  H2,  CO,  etc.,  and 

l 

j  =  WGS,  STR,  etc.  (82) 

The  spatial  distribution  of  oxygen  and  nitrogen  in  the  cathode  chan¬ 
nel  is  computed  analogously.  However,  the  mass  balance  equation 
for  oxygen  is  different  for  the  two  modeled  flow  patterns  for  the 
planar  cell  design,  namely  co-  and  counter-flow.  In  the  co-flow  case 
the  oxygen  content  decreases  with  the  axial  coordinate  when  a  cur¬ 
rent  is  produced.  For  the  counter-flow  case,  the  opposite  applies. 
The  different  boundary  conditions  for  the  considered  flow  patterns 
will  be  addressed  in  subsequent  chapters. 

dn^2'ca  -  0.0  (83) 

dx 


(^-  =  Rq2  with  Ro2  =  -ro2  ion  (co-flow)  or 
Ro2  =  ro2  ion  (counter-flow)  (84) 


Based  on  the  calculated  reaction  rates,  the  mass  balance  model  also 
computes  the  related  heat  source  term  according  to 

AHr  =  rHl  oxi  AHR)h2  oxi+rCOoxi  AHR  GOoxi+rdl-CH4  STR  AHR?  CH4  STR 

+rWGS  AHr?wgs  +  tREv— wgs  AHr?REV-wgs  (85) 

By  definition,  the  heats  of  reaction  are  attributed  to  the  solid  struc¬ 
ture  of  the  fuel  cell.  The  differing  heat  capacities  of  educt  and 
product  species  as  well  as  the  different  temperature  of  the  gas 
phase  and  the  solid  structure  result  in  an  enthalpy  flux  coupled 
to  the  mass  flow  from  the  solid  structure  to  the  gas  phase  and  vice 
versa.  Eqs.  (86)-(88)  give  the  mass  transfer  coupled  enthalpy  flux 
of  the  educts  and  products  at  the  anode  and  the  enthalpy  flux  at 
the  cathode: 


QsH,ed,an  —  ^  ^i,ed^p,i,TaK 

j  i 

(86) 

QsH,prod,an  —  ^  ^H,prod^p,i,TsK 

j  i 

(87) 

QsH,ed,ca  =  r02  ion^cI<Cp,02,TcK 

(88) 

In  these  equations,  r}  stands  for  the  reaction  rate  of  the  different 
considered  heterogeneous  reactions,  including  the  electrochemical 
reactions,  occurring  at  the  anode  and  r02  ion  represents  the  reaction 
rate  of  reaction  Eq.  (45). 

2.6.  Energy  balance  model 

Considering  the  strong  heat  exchange  mechanisms  between  the 
solid  cell  components,  i.e.  heat  conduction  and  radiation,  and  taking 
into  account  the  compactness  of  the  considered  fuel  cell  designs, 
the  assumption  of  only  one  effective  solid  structure  temperature  for 
all  solid  cell  components  is  reasonable  [11  ].  However,  this  assump¬ 
tion  does  not  necessarily  include  the  gas  phase  temperatures  due 
to  the  weak  heat  exchange  via  convection  and  the  negligible  radia¬ 
tive  heat  exchange  between  solid  structure  and  gas  phase.  In  order 
to  account  for  these  conclusions,  the  energy  balance  model  deter¬ 
mines  the  axial  temperature  profile  of  the  solid  structure,  of  the 
anode  gas  and  cathode  gas  as  well  as  the  air  in  the  air  delivery  tube 
in  case  of  a  tubular  design. 


Fig.  9.  Outline  of  the  energy  balance  model. 
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Fig.  9  depicts  the  outline  of  the  energy  balance  model  for  an 
infinitesimal  control  volume  with  the  three  balance  cases  of  the 
planar  design,  namely  the  anode  and  cathode  gas  channel  as  well 
as  the  solid  structure,  and  the  additional  fourth  balance  case  of 
tubular  cell  designs,  namely  the  air  delivery  tube. 

The  energy  balance  of  the  anode  gas  channel  includes  the  sen¬ 
sible  heat  flow,  han*Cp, an*Pai<»  entering  the  control  volume  at  the 
coordinate  x  and  leaving  it  at  x  +  dx.  Further,  the  enthalpy  flux  due  to 
the  different  heat  capacity  of  the  educt  and  product  species  reacting 
at  the  anode  are  accounted  for.  The  conductive  heat  stream  between 
solid  structure  and  anode  gas  channel  is  computed  via  Eq.  (89): 

Qconv.an  =  ^an^circ,an(^al<  —  7sl<)  with  OCa n  =  — -  (89) 

fthyd.an 

The  differential  equation  describing  the  axial  temperature  distri¬ 
bution  along  the  anode  gas  channel  can  thus  be  formulated  as 
follows: 

^'an(^an  al<^  =  OSH, prod, an  “  QsH.ed.an  “  Qconv.an  (90) 

It  should  be  noted,  that  this  model  considers  the  change  of  the 
total  molar  flow,  nan,  and  the  change  of  the  heat  capacity,  cP)an,  due 
to  the  change  of  the  gas  composition  over  the  control  volume.  In 
some  models  found  in  the  literature,  only  the  temperature  change 
is  considered  by  assuming  the  heat  capacity  and  sometimes  also 
the  molar  flow  are  fixed  at  the  inlet  conditions  of  the  according 
control  volume  [22,31,32].  For  coarse  discretization  aiming  at  short 
calculation  times,  this  can  affect  accuracy. 

The  energy  balance  of  the  solid  structure  involves  the  overall 
heat  of  reaction,  Eq.  (85),  and  the  mass  transfer  coupled  enthalpy 
fluxes  at  the  anode  and  cathode  electrode,  Eqs.  (86)-(88).  Besides, 
the  convective  heat  flows  between  the  solid  structure  and  anode 
as  well  as  cathode  channel  have  to  be  considered.  The  produced 
electrical  power  is  considered  as  source  term  for  the  solid  structure 
and  is  computed  according  to 


Pel  —  ~  4lact^tot^op  (91) 

The  convective  heat  flow  between  the  cathode  gas  channel  and  the 
solid  structure  is  determined  according  to  Eq.  (92).  The  sign  of  the 
convective  heat  flow  is  considered  positive,  when  heat  is  trans¬ 
ferred  from  the  solid  structure  to  the  cathode  gas.  This  reflects  the 
cooling  task  of  the  air  flow  in  the  cathode  channel: 

Qconv.ca  =  ^ca^circ.caC^sK  —  Tcl<)  with  OiCa  =  -  (92) 

fthyd.ca 

Finally,  the  energy  balance  of  the  solid  structure  includes  solid  heat 
conduction  yielding  a  second  order  derivate  differential  equation: 

d^T  i<  • 

^-s^cross  ^2  =  Qconv.ca  —  Qconv.an  —  QsH.ed.an  —  QsH.ed.ca 

+0sH,prod,an  +  AHr  +  Pej  (93) 

The  differential  equation  to  compute  the  cathode  gas  temperature 
of  the  planar  cell  designs  is  given  by 

d(hcaCp,ca7d<)  A  A  ,  A  ri 

- -  —  Qconv.ca  —  QsH.ed.an  +  Oconv,ADT  (CO-flow) 

=  QsH.ed.an  “  Qconv.ca  (counter-flow)  (94) 

In  principle,  the  tubular  cell  designs  have  a  co-flow  pattern.  EIow- 
ever,  the  design  typical  air  delivery  tube,  adds  one  more  term  to  the 
cathode  gas  energy  balance:  the  convective  heat  flow  from  the  air 
delivery  tube  to  the  cathode  gas,  0cOnv,ADT- 

Oconv.ADT  aggregates  the  convective  heat  transport  from  the 
cathode  air  to  the  air  delivery  tube  (ADT)  and  from  the  ADT  to 
the  cold  air  flowing  inside  the  ADT,  Eq.  (95).  In  this  model,  the 


Table  4 

Mass  balance  boundary  conditions 


Mass  balance 

Planar  co-flow 

Planar  counter-flow 

Tubular 

Anode  gas 

Cathode  gas 

Air  delivery  tube 

^an(O)  =  ^an,  0 
^ca(O)  =  ^ca,  0 

^an(O)  =  ^an,  0 
^ca(O)  =  ^ca,  0 

^an(O)  =  ^an,  0 
^catO)  =  'ica,  0 
^ADT(n)  =  ^ca,  0 

two  processes  are  accounted  for  via  an  effective  heat  exchange 
coefficient,  aADTieff,  which  was  proposed  in  [33].  The  calculation 
of  c^ADT  eff  includes  the  solid  heat  conduction  coefficient  of  the  ADT 
material  which  leads  to  either  enhanced  or  reduced  convective  heat 
exchange  in  cold  or  hot  regions  of  the  ADT,  respectively. 

Oconv.ADT  =  ^ADT,eff^circ,ADl(^adtK  “  ^cl<)  with  aADT,eff 


— ( - - - +—i —  In  ( r~P^)  + - - - \ 

ri,ADT  yH.ADT^ADT  4AS?ADT  y  H.ADT  J  ro,ADT^ca  J 


and  aADT  = 


Nu  Aadt 

dhyd,ADT 


(95) 


The  energy  balance  of  the  air  flow  inside  the  ADT  can  then  be  for¬ 
mulated  as  Eq.  (96).  Since  the  air  flow  inside  the  ADT  is  not  taking 
part  in  any  chemical  or  electrochemical  reactions,  the  molar  flow  is 
constant.  Accordingly  the  differential  equation  for  the  ADT  energy 
balance  was  simplified  in  the  present  model: 


d(cp,ADTradtI<)  . 

ft  ADT - -  =  Oconv.ADT 


(96) 


2.7.  Boundary  conditions 


The  mass  balance  boundary  condition  values  are  calculated 
based  on  the  user-defined  operating  conditions.  The  fuel  inlet  com¬ 
position  is  specified  in  molar  fractions  of  the  considered  species. 
The  total  molar  flow  of  the  fuel  gas  at  the  cell  inlet  is  either  calcu¬ 
lated  based  on  a  targeted  overall  current  density,  Jtot,av,  and  related 
fuel  utilization,  UF,  Eq.  (97),  or  by  specifying  a  lower  heating  value 
(LEIV)-based  input  power,  Pin: 

ftan.o  =  2F°UFPqH  with  EClH2  =yHl  +yco  +  4^ch4  (97) 


Tan.O 


Pin 

LHVinnRE 


(98) 


The  inlet  molar  flow  of  the  cathode  gas  is  either  assigned  directly 
or  by  specifying  an  air-to-fuel  ratio,  X: 


ftca.O 


ft02,stoic^ 

y°2 


with  n02, stoic 


=  han,o(0-5yH2  +  0.5y<;o  +  2ycH4)  (99) 

Depending  on  the  investigated  fuel  cell  design  and  the  considered 
flow  pattern,  the  mass  balance  boundary  values  were  assigned  to 
different  control  volumes  (CV),  see  Table  4  where  (0)  denotes  the 
first  and  (n)  the  last  CV.  The  total  molar  flow  values  nan,  ftca  and 
nADT  are  related  to  the  molar  flows  of  the  single  species  via  the 
user-defined  molar  fractions. 

Besides  the  inlet  molar  flow,  the  present  model  requires  the  def¬ 
inition  of  the  gas  temperatures  at  the  inlet  in  order  to  solve  the 


Table  5 

Gas  phase  energy  balance  boundary  conditions 


Energy  balance 

Planar  co-flow 

Planar  counter-flow 

Tubular 

Anode  gas 

^aK(O)  =  Tm,in 

Tal<(0)  =  Tan, in 

Tal<(0)  =  Tan, in 

Cathode  gas 

T:I<(0)  =  Tea, in 

Tel<(n)  =  Tea, in 

T:I<(0)  =  Tadtl<(0) 

Air  delivery  tube 

- 

- 

Tadtl<(n)  =  ^ADT.in 

Table  6 

Fuel  and  cathode  gas  compositions  of  the  BMT 
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Table  7 

Geometrical  input  data  of  BMT 


energy  balance  equation  of  the  gas  channels.  The  inlet  gas  temper¬ 
atures  of  the  anode,  the  cathode  and  eventually  of  the  air  deliver 
tube  are  user-defined  values  and  assigned  as  shown  in  Table  5.  In 
case  of  the  tubular  cell  design,  the  boundary  conditions  of  the  cath¬ 
ode  gas  channel  ensure  the  continuity  between  the  temperature 
and  the  molar  flow  of  the  air  flow  exiting  the  air  delivery  tube  and 
entering  the  cathode  gas  channel. 

Owing  to  the  considered  solid  heat  conduction,  the  energy  bal¬ 
ance  of  the  solid  structure  includes  a  second  derivate  term.  Hence, 
the  solution  of  each  CV  depends  on  the  solutions  of  the  neighboring 
CVs.  Consequently,  boundary  conditions  have  to  be  defined  at  both 
ends  of  the  integration  region.  The  differential  equation  applying 
at  the  left  boundary  (first  CV  of  the  integration  region)  is  given  by 

AsACross^J=Qhloss  (100) 

Qhioss  represents  the  heat  loss  flow  to  the  surroundings.  At  the  right 
boundary,  the  source  terms  of  the  last  CV  of  the  integration  region 
have  to  be  taken  into  account: 

Across  ^  =  (Qconv.ca  —  Qconv.an  —  0sH,ed,an  —  0sH,ed,ca 

+0sH,prod,an  +  AHr  +Pei)dx  —  Qhioss  (101 ) 

For  the  tubular  design  the  cell  ends  are  generally  assumed  adiabatic. 
In  contrast,  our  model  features  the  possibility  to  assume  adiabatic 
cell  ends  or  conductive  heat  losses  at  the  cell  ends  of  planar  cells. 
Where  adiabatic  cell  ends  are  assumed  for  the  simulation  of  a  per¬ 
fect  insulation,  Qhi0Ss  is  set  to  zero.  Assuming  conductive  heat  losses 
through  insulation,  Eq.  (102)  applies: 

Qhioss  —  ^insulAt/p,  cross  (^sK  —  lamb)  (102) 


3.  Model  validation 

3.1  Planar  geometry 

A  verification  of  our  model  is  carried  out  by  comparing  the  com¬ 
puted  results  to  a  benchmark  test  (BMT).  The  BMT  was  defined 
in  an  International  Energy  Agency  (IEA)  program  for  the  numer¬ 
ical  simulation  of  SOFCs  with  planar  geometry  [34].  In  total,  nine 
different  academic  and  industrial  institutions  participated  in  the 
BMT.  The  independently  developed  models  showed  good  agree¬ 
ment  concerning  the  predicted  physical  behavior  for  two  fuel  gas 
compositions,  in  detail  humidified  hydrogen  and  30%  pre-reformed 
methane,  Tables  6  and  7  gives  the  BMT  macro  and  micro  geometry 
data  of  the  considered  planar,  electrolyte  supported  SOFC.  Table  8 
shows  the  operational  conditions  of  the  considered  SOFC  defined 
by  the  BMT  participants. 

The  BMT  participants  agreed  to  consider  only  ohmic  losses  for  a 
better  comparability  of  the  model  results.  The  specific  ohmic  resis¬ 
tances  of  the  ceramic  components  were  calculated  according  to  Eqs. 
(4)-(7).  For  the  planar  SOFC,  ceramic  bipolar  plates  were  assumed. 
Further,  the  applied  kinetics  data  for  the  steam  reforming  reaction 


Table  8 

Operational  conditions  of  BMT 


Model  input  data 

Unit 

Value 

Source 

Operational  conditions 

Targeted  mean  current  density 

Am-2 

3000.0 

Targeted  fuel  utilization 

% 

85.0 

Air-to-fuel  ratio 

- 

7.0 

Fuel  gas  inlet  temperature 

I< 

1173.15 

[34] 

Cathode  gas  inlet  temperature 

I< 

1173.15 

Ambient  temperature 

I< 

293.15 

System  pressure 

bara 

1.01325 
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Table  9 

Energy  balance  input  data  of  BMT 


Model  input  data 

Unit 

Value 

Source 

Solid  heat  conduction 

Heat  conductivity  of  anode 

2.0 

Heat  conductivity  of  electrolyte 

Wm1  K-1 

2.0 

[34] 

Heat  conductivity  of  cathode 

2.0 

Heat  conductivity  of  interconnect 

2.0 

Convective  heat  transfer 

Nusselt  number 

- 

4.0 

[34] 

(STR)  were  prescribed: 

rcH4STR,Ach  =  4274.0  mol  m“2  s_1  bara-1  pCH4 

(  -82000J  mol-1  \  Mno, 

xexp(, — — )  (,03) 

Table  9  gives  the  BMT  parameter  used  for  the  energy  balance 
calculations.  Besides,  several  assumptions  were  defined  for  the 
BMT,  Table  10.  The  most  important  assumptions  were  that  the 
activation  polarization  losses  are  equal  to  the  ohmic  loss  of  the 
electrolyte  and  that  diffusion  losses  are  negligible.  For  this  rea¬ 
son,  a  validation  of  activation  and  diffusion  loss  model  parameters 
is  not  possible  using  the  BMT.  Furthermore,  the  change  of  the 
molar  flow  and  the  heat  capacity  of  the  gas  mixture  was  neglected 
in  the  calculation  of  the  gas  phase  energy  balance.  Instead  the 
calculations  were  performed  with  a  simplified  definition  of  the 


Table  10 

Model  settings  of  BMT 


Model  settings 

Description 

Electrochemical  loss  model  settings 

Considered  electrochemical  active  species 

Hydrogen 

Activation  polarization  equation 

Equal  to  ohmic  loss  of 
electrolyte 

Energy  balance  settings 

Solid  heat  transfer  mechanism 

Non-isothermal  with  solid 
heat  conduction 

Coupled  heat  and  mass  transport 

Considered 

Heat  loss  mechanism  at  outer  surface 

Conduction  through 
insulating  plates 

Heat  capacity  and  heat  of  reaction  correlation 

Published  in  [11] 

Definition  of  sensible  heat  gradient 

ncp|x.f 

sensible  heat  gradient  by  holding  the  molar  flow  and  heat  capac¬ 
ity  constant  at  the  value  of  the  antecedent  control  volume  for  the 
respective  control  volume.  Further  information  can  be  found  in 
[34]. 

Fig.  10  shows  the  comparison  of  the  results  of  our  model  for  the 
co-flow  pattern  considering  both  IEA  fuel  gases,  Table  6.  The  terms 
“Minimum  value”  and  “Maximum  value”  refer  to  the  minimum  and 
maximum  values  among  all  the  BMT  participants,  whereas  “Own” 
refers  to  the  results  obtained  from  our  model.  Fig.  11  shows  the 
comparison  of  the  results  for  the  counter-flow  pattern.  It  can  be 
seen,  that  our  model  behaves  physically  correct  and  that  there  is 
no  systematic  discrepancy  to  the  results  of  other  participants  in 


Planar  co-flow  model  comparison  for  IEA  benchmark  1  input 


4500 
w  4000 
§  3500 
?  3000 
?  2500 

Q.  2000 

(/> 

£  1500 
o  1000 
500 
0 


□  Minimum  value  □  Own  ■  Maximun  value 


Planar  co-flow  model  comparison  for  IEA  benchmark  2  input 
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Fig.  10.  Comparison  of  co-flow  model  results  with  BMT  results  for  fuel  compositions  IEA  1  and  IEA  2. 
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Planar  counter-flow  model  comparison  for  IEA  benchmark  1  input 
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Planar  counter-flow  model  comparison  for  IEA  benchmark  2  input 
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Fig.  11.  Comparison  of  counter-flow  model  results  with  BMT  results  for  fuel  compositions  IEA  1  and  IEA  2. 


the  BMT.  It  is  concluded,  that  our  model  is  a  reliable  tool  for  the 
simulation  of  planar  SOFCs. 

3.2.  Tubular  geometry 

The  aim  of  validating  our  model  for  the  standard  tubular  cell 
geometry  is  to  assure  that  the  generalized  approach,  which  is  based 
on  characteristic  lengths  and  areas,  has  been  properly  adapted  to 
the  tubular  cell  design.  Although  the  most  important  parts  of  the 
model  code  were  already  verified  using  the  BMT  for  the  planar 
geometry,  there  is  still  a  need  to  also  check  the  equations  involved  in 
the  calculation  of  the  activation  and  diffusion  losses  as  those  were 
excluded  from  the  BMT.  The  geometric  data  for  the  description  of 
the  standard  tubular  cell  is  given  in  Table  11. 

The  experimental  voltage-current  curves  used  for  the  valida¬ 
tion  were  published  in  [35].  The  characterized  standard  tubular 
cell  was  operated  with  a  fuel  gas  mixture  consisting  of  89vol.% 
H2  and  11  vol.%  H20  at  three  different  cell  temperatures,  900  °C, 
940  °C  and  1000  °C  with  a  fuel  utilization  of  85%  and  an  air-to- 
fuel  ratio  of  4.  The  latter  feature  allows  checking  for  the  correct 
prediction  of  temperature  effects  and  hence  makes  this  data  set 
predestined  for  model  validation  purpose.  Model  results  obtained 
for  other  temperatures  can  be  assumed  valid,  in  case  the  model  reli¬ 
ably  predicts  the  temperature  induced  effects  of  the  validation  case. 
When  model  validation  is  carried  out  with  data  sets  measured  at  a 
single  temperature,  the  results  are  only  valid  for  this  temperature 
and  extrapolated  results  have  to  be  considered  less  reliable. 


For  the  simulation,  the  gas  inlet  temperatures  were  assumed 
100 K  lower  than  the  respective  mean  cell  temperature  [28].  The 
heat  transfer  calculation  input  data  is  given  in  Table  12.  The  elec¬ 
trochemical  loss  model  parameters  are  given  in  Table  13. 

The  conductivity  of  the  anode,  electrolyte  and  cathode  are  cal¬ 
culated  using  the  same  correlations  as  in  the  BMT.  The  diffusion 
polarization  loss  parameters  were  taken  from  Stiller  [31].  FIow- 
ever,  the  cathode  electrode  tortuosity  was  slightly  increased  from 


Table  11 

Geometrical  input  data  for  standard  tubular  cell  model 


Model  input  data 

Unit 

Value 

Source 

Macrogeometry  of  tubular  cell 

Flow  design 

Co-flow  with 
air  delivery 
tube 

Cell  tube  length 

m 

1.5 

[28,31] 

Inner  radius  of  air  delivery  tube 

m 

0.0025 

Outer  radius  of  air  delivery  tube 

m 

0.004 

Inner  radius  of  cell  tube 

m 

0.00866 

Thickness  of  virtual  fuel  channel 

m 

0.0023 

(depends  on  stacking  of  cell  tubes) 

Percentage  of  circumferential  length 

% 

10.0 

of  cell  tube  covered  by  IC 

Microgeometry  of  planar  cell 

Support  design 

- 

Cathode 

Anode  thickness 

pan 

100.0 

[28,31] 

Electrolyte  thickness 

pan 

40.0 

Cathode  thickness 

pan 

2200.0 
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Table  12 

Energy  balance  input  data  for  standard  tubular  cell  model  validation 


Model  input  data 

Unit  Value 

Source 

Solid  heat  conduction 

Heat  conductivity  of  anode 

6.23 

Heat  conductivity  of  electrolyte 

i  1  2.7 

Heat  conductivity  of  cathode 

Wm-1K-1  g6 

[31] 

Heat  conductivity  of  air  delivery  tube 

11.8 

Nusselt  number 

4.0 

Table  13 

Electrochemical  loss  model  input  data  for  tubular  and  Delta8  cell  model 


Model  input  data 

Unit 

Value 

Source 

Activation  polarization 

H2  oxidation  activation  energy 

J  mol-1 

110000.0 

[13] 

H2  oxidation  pre-exponential  factor 

A  ITT2 

7000000000.0 

CO  oxidation  activation  energy 

J  mol-1 

110000.0 

[11,13] 

CO  oxidation  pre-exponential  factor 

Am-2 

5000000000.0 

02  reduction  activation  energy 

J  mol-1 

149500.0 

[28] 

02  reduction  pre-exponential  factor 

Am-2 

10260000000.0 

Ohmic  polarization 

Electric  conductivity  of  anode 

1  m-1 

T-dependent, 

Eq.  (4) 

[12] 

Ionic  conductivity  of  electrolyte 

1  m-1 

T-dependent, 

Eq.  (5) 

Electric  conductivity  of  cathode 

1  fi-1  m-1 

T-dependent, 

Eq.  (6) 

Electric  conductivity  of  interconnect 

Q,  m 

2.0E-7 

[28] 

Diffusion  polarization 

Porosity  of  anode 

- 

0.4 

Tortuosity  of  anode 

- 

3.0 

[31] 

Pore  diameter  of  anode  and  cathode 

m 

1.0E-6 

Porosity  of  cathode 

- 

0.5 

Tortuosity  of  cathode 

- 

1.8 

- 

1.5  to  1.8  in  order  to  better  reproduce  the  diffusion  limitation  at 
high  current  densities.  The  activation  polarization  parameters  for 
the  hydrogen  oxidation  and  for  the  oxygen  ionization  were  taken 
from  Campanari  [13],  and  Thorud  [28],  respectively.  The  activation 
loss  parameters  for  the  carbon  monoxide  oxidation  were  calculated 
based  on  those  of  the  hydrogen  oxidation  and  the  relation  proposed 
in  [11]. 

In  the  following  the  model  results  are  compared  to  the  experi¬ 
mental  voltage-current  curves.  As  the  same  voltage-current  curves 
were  used  for  model  validation  purpose  in  [28],  the  detailed  discus¬ 
sion  of  the  reasons  for  the  difference  between  the  model  results  and 
the  measurement,  carried  out  in  [28],  is  not  repeated  at  this  point. 
However,  the  higher  accuracy  of  the  present  model  is  highlighted. 


Current  density  [A/m2] 

— ■— 1000°C  Measurement  — •—  940°C  Measurement  90CTC  Measurement 

1000*0  Simulation  — 940*0  Simulation _ 900°C  Simulation 


Fig.  12.  Comparison  of  simulated  and  measured  voltage-current  curve  of  a  standard 
tubular  cell  measured  at  mean  cell  temperature  of  900  °C,  940  °C  and  1000  °C. 


Fig.  12  shows  the  comparison  of  simulated  (dashed  line)  and 
measured  (full  line)  voltage-current  density  (U-I)  characteristics 
of  the  tubular  cell  operated  at  900  °C,  940 °C  and  1000  °C.  The 
mean  variance  of  the  predicted  voltage  for  900  °C  operational  mean 
cell  temperature  compared  to  the  measured  values  is  1.6%  and  the 
maximum  variance  is  below  ±3%.  For  940  °C  the  predicted  voltage 
values  diverge  by  less  than  1  %  for  small  and  medium  current  density 
values.  At  higher  current  density  values  the  variance  increases  up 
to  5%.  The  mean  variance  between  predicted  and  measured  volt¬ 
age  values  for  940  °C  mean  cell  temperature  is  1.3%.  At  a  mean 
cell  temperature  of  1000  °C,  the  mean  variance  between  predicted 
and  measured  voltage  values  is  2.3%  and  the  maximum  variance 
is  below  ±4%.  In  contrast  to  [28],  at  high  current  density  values, 
the  present  model  predicts  a  slight  bend  of  the  U-I  curve  due 
to  increasing  diffusion  losses.  This  phenomenon  can  be  observed 
more  pronounced  for  the  measurement.  The  prediction  of  this  phe¬ 
nomenon  emphasizes  that  our  tubular  model  behaves  physically 
correct. 

It  is  concluded,  that  our  generalized  model  was  properly  adapted 
to  the  tubular  cell  geometry.  Further,  the  activation  and  diffusion 
loss  model  equations  and  related  model  parameters  proved  to  yield 
accurate  results. 

3.3.  Triangular  geometry 

The  published  data  dealing  with  Delta8  cells  (D8)  is  scarce  as 
they  have  only  recently  been  developed.  Nevertheless,  Table  14 
gives  the  available  D8  cell  geometrical  data.  It  is  important  to  point 
out  that  the  geometry  of  the  air  delivery  tube  was  not  taken  from 
a  publication  but  was  assumed  to  yield  an  adequate  cathode  gas 
channel  cross-sectional  area.  This  assumption  directly  impacts  the 
convective  heat  transfer  between  the  cathode  gas  and  the  air  deliv¬ 
ery  tube.  The  convective  heat  transfer  between  the  cathode  gas  and 
the  solid  structure  is  however  dominating  in  any  case,  such  that  the 
overall  impact  of  this  assumption  is  considered  less  significant. 

The  experimental  voltage-current  curve  used  for  the  validation 
of  the  triangular  cell  model  was  measured  using  a  D8  cell  segment 
[3].  The  1  cm  long  D8  cell  segment  was  operated  under  isother¬ 
mal  conditions  in  an  oven  at  temperatures  between  950  °C  and 
965  °C  with  fuel  and  air  excess.  The  fuel  mixture  consisted  of  50 
vol.%  H2  and  50vol.%  H20,  being  the  expected  average  gas  com¬ 
position  of  a  full-scale  cell  operated  at  85%  fuel  utilization  with 
89 vol.%  H2  and  11  vol.%  H20.  Due  to  the  isothermal  conditions 
and  the  virtually  constant  reactant  partial  pressures  along  the  D8 
cell  segment,  the  voltage-current  curve  was  simulated  by  using 
the  electrochemical  performance  model  without  considering  the 

Table  14 

Geometrical  input  data  for  Delta8  cell  model 
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Fig.  13.  Comparison  of  simulated  and  measured  voltage-current  curve  of  a  D8  cell 
segment  measured  at  a  temperature  between  950  °C  and  965  °C. 

energy  and  heat  balance.  Concerning  the  model  parameters,  it  was 
assumed  that  the  anode  and  cathode  materials  used  for  the  D8 
cells  are  like  those  used  for  the  standard  tubular  cells.  Accordingly, 
the  model  parameters  given  in  Table  13  apply.  Exceptions  were 
the  electric  conductivity  of  the  interconnect  material,  which  was 
calculated  according  to  Eq.  (7),  and  the  ionic  conductivity  of  the 
electrolyte,  as  D8  cells  employ  scandia-stabilized  zirconia  (ScSZ) 
instead  of  the  state-of-the-art  yttria-stabilized  zirconia  (YSZ)  as 
electrolyte  material  [3].  At  temperatures  above  900  °C,  ScSZ  fea¬ 
tures  ionic  conductivities  about  three  times  higher  than  YSZ  [36]. 
Hence,  the  conductivity  values  computed  via  Eq.  (5)  were  multi¬ 
plied  by  three  for  the  D8  cell  segment  simulation.  Fig.  13  shows  the 
comparison  of  the  measured  and  the  simulated  voltage-current 
curve  of  the  D8  cell  segment.  The  simulation  result  for  957  °C  fits 
the  measured  data  very  well.  It  is  concluded,  that  the  present  model 
was  successfully  applied  to  the  triangular  geometry  and  that  the 
above  discussed  assumptions  are  valid. 

4.  Performance  comparison 

Benchmarking  of  SOFCs  is  usually  performed  with  humidified 
hydrogen.  The  fuel  of  choice  for  SOFCs  is  however  pre-reformed  nat¬ 
ural  gas,  which,  besides  electrochemical  reactions,  involves  steam 
reforming  reactions  and  enhanced  heat  transfer  processes.  In  order 
to  account  for  the  whole  variety  of  possible  interactions  between 
charge,  heat  and  mass  transfer  processes,  the  comparison  of  the 


Table  15 

Electrochemical  loss  model  input  data  for  planar  cell  model 


Model  input  data 

Unit 

Value 

Source 

Activation  polarization 

H2  oxidation  activation  energy 

J  mol-1 

120000.0 

[22] 

H2  oxidation  pre-exponential  factor 

Am-2 

2900000000.0 

CO  oxidation  activation  energy 

J  mol-1 

120000.0 

[11] 

CO  oxidation  pre-exponential  factor 

Am-2 

2070000000.0 

02  reduction  activation  energy 

J  mol-1 

120000.0 

[22] 

02  reduction  pre-exponential  factor 

Am-2 

7000000000.0 

Ohmic  polarization 

Electric  conductivity  of  anode 

1  ft”1  m-1 

T-dependent,  Eq.  (4) 

Ionic  conductivity  of  electrolyte 

1  ft”1  m-1 

T-dependent,  Eq.  (5) 

[12] 

Electric  conductivity  of  cathode 

1  ft-1  m_1 

T-dependent,  Eq.  (6) 

Electric  conductivity  of  interconnect 

1  ft-1  m_1 

T-dependent,  Eq.  (7) 

Diffusion  polarization 

Porosity  of  anode 

- 

0.5 

Tortuosity  of  anode 

- 

3.0 

Pore  diameter  of  anode  and  cathode 

m 

1.0E-6 

[10] 

Porosity  of  cathode 

- 

0.5 

Tortuosity  of  cathode 

- 

3.0 

Cell  length  [-] 

- Planar  •  -  -  -  Tubular - Delta8| 


Fig.  14.  Current  density  distribution  of  planar,  tubular  and  Delta8  cell. 

planar,  tubular  and  Delta8  design  was  conducted  considering  pre¬ 
reformed  NG  as  fuel  gas  The  electrochemical  performance  model 
input  parameters  employed  in  the  planar  model  are  summarized 
in  Table  15.  Note  that,  as  only  the  conversion  of  hydrogen  was 
considered  in  [22],  no  parameters  for  the  carbon  monoxide  acti¬ 
vation  polarization  were  specified.  However,  in  [11]  it  is  stated 
that  the  activation  potential  loss  of  the  electrochemical  conversion 
of  carbon  monoxide  is  approximately  1.4  times  higher  than  that 
for  the  electrochemical  hydrogen  conversion.  Based  on  this,  the 
pre-exponential  factor  for  the  calculation  of  the  carbon  monoxide 
conversion  exchange  current  density  was  determined. 

For  a  better  comparability  of  the  partial  pressure  related  voltage 
losses,  the  fuel  utilization  and  the  operational  voltage  were  fixed  at 
85%  and  0.6  V,  respectively.  Further,  the  air-to-fuel  ratio  was  fixed 
at  a  value  of  4.28.  The  fuel  and  air  inlet  temperatures  were  adjusted 
to  meet  the  typical  mean  cell  temperature  of  the  investigated  cell 
designs  being  950  °C  for  the  planar  and  Delta8  cell  and  1000  °C  for 
the  tubular  cell. 

Fig.  14  shows  the  current  density  plotted  against  the  dimension¬ 
less  cell  length  of  the  three  investigated  cell  designs.  In  average,  the 
electrolyte  supported  planar  cell  produces  3468  Am-2  while  the 
tubular  cell  yields  2890  Am-2.  The  Delta8  cell  has  the  highest  cur¬ 
rent  output  with  5027  Am-2,  which  corresponds  to  800  W  at  0.6  V 
operational  voltage.  This  value  is  close  to  the  projected  power  out¬ 
put  of  1000  W  at  a  mean  cell  temperature  of  1000  °C  [3],  however 
some  improvements  will  still  be  necessary  to  reach  this  target  in 
case  it  was  defined  for  natural  gas  operation  of  the  Delta8  cells. 

The  shape  of  the  current  density  distribution  is  similar  for 
the  tubular  and  the  Delta8  cell.  The  maximum  current  density  is 
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Fig.  15.  Solid  and  cathode  air  temperature  distribution  of  planar,  tubular  and  Delta8 
cell. 
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Fig.  16.  Methane  molar  fraction  distribution  of  planar,  tubular  and  Delta8  cell. 

reached  in  the  first  half  of  the  cell  length,  while  the  current  density 
maximum  of  the  planar  cell  can  be  observed  in  the  second  half. 
The  reason  for  this  can  be  found  in  the  temperature  distribution. 
Fig.  15  depicts  the  temperature  distribution  over  the  dimension¬ 
less  cell  length.  In  contrast  to  the  tubular  and  the  Delta8  cell,  which 
feature  an  air  delivery  tube  cooling  the  cell  outlet  region,  the  solid 
temperature  of  the  planar  cell  keeps  increasing  towards  the  cell  end. 

Due  to  the  short  cell  length  of  the  planar  cell,  the  planar  cell 
design  features  the  highest  thermal  gradients  with  maximum  val¬ 
ues  as  high  as  3.3  K  mm-1.  The  tubular  cell  exhibits  lower  thermal 
gradients  of  at  most  1.3  K  mm-1.  Despite  the  fact  that  the  Delta8 
cell  was  assumed  to  have  a  lower  solid  heat  conductivity  coefficient 
of  2  W  m-1  K-1  compared  to  9.3  W  m-1  K-1  of  the  tubular  cell,  the 
Delta8  cell  shows  the  lowest  thermal  gradients  with  0.8  W  m-1  K-1 . 
This  can  be  assigned  to  the  enhanced  convective  heat  exchange 
between  the  solid  structure  and  the  cathode  air  due  to  a  smaller 
hydraulic  diameter  and  the  relatively  seen  bigger  heat  exchange 
surface  of  the  cathode  gas  channel.  The  same  reason  applies  for 
the  less  pronounced  reforming  cold  spot  of  the  Delta8  cell,  which  is 
inhibited  by  the  pre-heated  cathode  air  exiting  the  air  delivery  tube 
providing  heat  for  the  endothermal  reforming  reactions.  This  heat 
exchange  between  the  cathode  air  and  the  solid  structure  of  the 
tubular  cell  is  limited  yielding  a  considerable  cold  spot  at  the  cell 
inlet  and  higher  thermal  gradients.  Another  reason  for  the  lower 
thermal  gradients  of  the  Delta8  cell  can  be  found  in  the  steam 
reforming  reactions. 

Fig.  16  shows  the  methane  molar  fraction  as  a  function  of  the 
absolute  cell  length.  The  endothermal  reforming  reactions  (STR) 
entail  almost  constant  solid  temperatures  for  all  the  investigated 
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Fig.  18.  Voltage  loss  distribution  of  tubular  cell. 

cell  types  until  approximately  70%  of  the  methane  is  converted. 
The  planar  cell  exhibits  an  approximately  six  times  lower  cata¬ 
lyst  load  and  the  highest  cell  inlet  temperature.  In  consequence, 
the  STR  is  complete  after  3  cm  of  the  cell  length.  The  tubular  and 
the  Delta8  cell  have  approximately  the  same  catalyst  load.  Despite 
about  50 1<  higher  inlet  temperatures,  the  STR  occurs  slower  in  the 
Delta8  cell  being  complete  after  approximately  36  cm  versus  30  cm 
for  the  tubular  cell.  The  reasons  for  that  are  the  more  pronounced 
diffusion  limitation  due  to  the  almost  doubled  hydraulic  diameter 
of  the  anode  channel  and  an  approximately  30%  higher  catalyst  load 
of  the  Delta8  cell  compared  to  the  tubular  cell. 

The  local  distribution  of  the  voltage  losses  are  depicted  in  Fig.  17 
for  the  planar  geometry,  in  Fig.  18  for  the  tubular  geometry  and  in 
Fig.  19  for  the  Delta8  geometry,  respectively.  Note  that  the  volt¬ 
age  losses  are  plotted  is  in  logarithmic  scale  on  the  ordinate.  For 
the  planar  geometry,  diffusion  losses  are  approximately  two  to 
three  orders  of  magnitude  lower  than  all  other  voltage  losses.  This 
is  due  to  the  very  short  diffusion  path  length  of  the  investigated 
electrolyte  supported  cell  micro  geometry.  Concerning  the  anode 
activation  losses,  the  planar  cell  differs  significantly  from  the  tubu¬ 
lar  and  Delta8  cell.  In  the  front  parts  of  the  planar  cell,  anode  and 
cathode  activation  losses  are  on  comparable  level  with  the  ohmic 
losses.  Towards  the  cell  end  the  cathode  activation  and  the  ohmic 
losses  decrease  strongly  yielding  values  about  one  order  of  magni¬ 
tude  smaller  than  the  anode  activation  losses.  This  can  be  explained 
by  the  high  prevalent  temperature  at  the  cell  outlet  which  induces 
high  conductivity  of  the  YSZ  electrolyte.  The  decreasing  anode 
product  partial  pressures  towards  the  cell  end,  due  to  fuel  deple- 
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Fig.  17.  Voltage  loss  distribution  of  planar  cell. 


Fig.  19.  Voltage  loss  distribution  of  Delta8  cell. 
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Fig.  20.  Area  specific  resistance  distribution  of  planar,  tubular  and  Delta8  cell. 

tion,  in  contrast  inhibit  such  a  decrease  for  the  anode  activation 
losses.  This  phenomenon  can  neither  be  observed  for  the  tubular 
nor  for  the  Delta8  cell,  indicating  that  the  employed  model  param¬ 
eters  characterize  an  electrochemically  very  active  anode  catalyst 
used  in  these  cells  which  allows  electrochemical  conversion  with 
very  low  voltage  losses  even  at  low  educt  partial  pressures.  As  a 
result  of  the  assumption  that  the  tubular  and  the  Delta8  cell  fea¬ 
ture  the  same  anode  and  cathode  electrode  material,  the  voltage 
loss  pattern  of  both  cells  is  quite  similar.  Anodic  diffusion  losses 
are  negligible  as  is  also  the  case  for  the  planar  cell  design.  In  the 
cell  regions  where  the  steam  reforming  reactions  cause  very  low 
solid  temperatures,  the  cathode  activation  losses  are  higher  than 
the  ohmic  losses.  The  picture  changes  however  after  completion 
of  the  reforming  reactions  with  the  associated  strong  temperature 
increase.  In  case  of  the  tubular  cell,  the  cathode  activation  losses 
decrease  to  values  even  lower  than  the  cathode  diffusion  losses.  A 
similar  trend  can  be  observed  for  the  Delta8  cell.  However,  the  thin¬ 
ner  cathode  electrode  of  the  Delta8  geometry  inhibits  the  diffusion 
losses  to  be  come  higher  than  the  cathode  activation  losses. 

The  ohmic  losses  are  very  important  for  all  three  investigated 
cell  designs.  Ohmic  losses  are  straightforwardly  computed  via  the 
prevalent  current  density  and  the  area  specific  resistance.  The  area 
specific  resistance  depends  on  the  temperature-dependent  ionic 
and  electronic  conductivity  and  the  geometry  of  the  current  con¬ 
ducting  parts  of  the  cells.  The  calculation  was  discussed  in  detail 
before.  Fig.  20  shows  the  area  specific  ohmic  resistance  (ASR)  of  the 
three  investigated  cells  plotted  vs.  the  dimensionless  cell  length. 
The  tubular  cell  features  the  highest  ASR  due  to  the  longest  cur¬ 
rent  path  with  approximately  31  mm.  Further,  the  tubular  cell  ASR 
is  strongly  linked  to  the  prevalent  solid  temperature  due  to  the 
employed  YSZ  electrolyte.  In  the  cold  cell  inlet  region  where  the 
cell  temperature  is  about  1000  K,  the  ionic  conductivity  of  the  YSZ 
electrolyte  is  only  one  tenth  of  value  corresponding  to  the  mean 
cell  temperature  of  1273  K.  The  strong  temperature  increase  after 
completion  of  the  steam  reforming  reactions  explains  the  steep 
decrease  of  the  ASR  of  the  tubular  cell.  The  Delta8  cell  in  contrast 
shows  only  a  slight  decrease  of  the  ASR  despite  that  the  local  cell 
temperatures  increase  from  1050 1<  at  the  cell  inlet  to  1300  K  at  the 
cell  outlet.  The  reason  for  this  is  that  the  dominating  role  of  the 
electrolyte  for  the  ASR  is  weakened  by  the  employed  ScSZ  elec¬ 
trolyte  material.  The  rest  of  the  employed  conducting  materials, 
namely  the  electrode  and  interconnect  materials,  are  not  as  sen¬ 
sitive  towards  the  local  temperature.  Despite  that  the  Delta8  cell 
features  a  ScSZ  electrolyte  with  approximately  one  third  of  the  spe¬ 
cific  resistance  of  the  YSZ  electrolyte  and  a  20%  shorter  current  path 
length,  the  corresponding  ASR  is  in  average  only  50%  lower  than  that 
of  the  tubular  cell.  The  reason  for  that  can  be  found  in  the  30%  thin¬ 


ner  cathode  electrode  of  the  Delta8  cell  and  the  according  reduction 
of  the  current  conducting  cross-sectional  area  which  compensates 
the  shorter  current  path. 

As  expected,  the  ASR  of  the  planar  cell  is  significantly  lower  than 
that  of  the  tubular  cell  due  to  its  very  short  current  path  length  of 
approximately  3  mm.  However,  the  thick  YSZ  electrolyte  induces  a 
strong  temperature  dependency  of  the  planar  cell  ASR  which  leads 
to  higher  values  than  the  Delta8  cell  at  the  cell  inlet  despite  approx¬ 
imately  50 1<  higher  solid  temperatures.  Comparing  the  Delta8  cell 
with  the  planar  cell  shows  that  in  average  the  ASR  of  both  cell 
designs  is  on  a  comparable  level.  That  clearly  demonstrates  the 
improvements  which  were  made  in  the  development  of  the  Delta8 
cell  compared  to  its  predecessor,  the  standard  tubular  cell  design. 

5.  Conclusions 

A  generalized,  finite  volume-based  SOFC  model  has  been  devel¬ 
oped.  The  model  includes  charge,  mass  and  heat  transport  as  well 
as  a  Langmuir-Hinshelwood  type  applied  kinetics  models  for  the 
steam  reforming  reaction.  The  model  development  was  focused 
on  a  fast  applicability  to  various  cell  geometries,  short  calcula¬ 
tion  times  to  allow  for  system  analysis  calculations  and  high  fuel 
flexibility.  The  model  was  applied  to  a  planar,  the  standard  tubu¬ 
lar  and  the  triangular  tube  cell  geometry  (Delta8),  showing  good 
agreement  with  experimental  and  benchmark  test  data.  Based  on 
the  computed  local  species,  temperature,  current  density  and  volt¬ 
age  loss  distributions,  the  performance  of  the  three  investigated 
cell  designs  was  discussed,  giving  deep  insight  into  the  interrela¬ 
tions  between  cell  geometry  and  transport  phenomena  inside  the 
cells.  It  was  shown  that  the  performance  of  the  tubular  design  was 
improved  by  introducing  the  Delta8  designs.  The  Delta8  design 
features  a  more  homogeneous  temperature  distribution  than  the 
standard  tubular  cell  due  to  improved  heat  exchange  between  cath¬ 
ode  air  and  solid  structure.  This  eliminates  the  strong  cold  spot  at 
the  cell  inlet  observed  for  the  standard  tubular  cell.  This,  together 
with  the  employed  ScSZ  electrolyte,  results  in  area  specific  resis¬ 
tance  values  of  the  Delta8  cell  which  are  approximately  half  of  the 
standard  tubular  values.  In  consequence,  the  performance  of  the 
new  Delta8  cell  design  is  almost  double  that  of  the  standard  tubular 
design.  The  planar  electrolyte  supported  cell  is  also  outperformed 
by  the  Delta8  cell.  In  future,  our  model  will  be  included  in  a  flow¬ 
sheeting  software  package  allowing  for  detailed  systems  analysis 
calculations  of  SOFC-based  power  plants  applying  different  fuels 
such  as,  e.g.  biomass  gasification  derived  producer  gases,  biogas, 
coal  gas,  etc. 
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